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ABSTRACT 


In  this  thesis  we  derive  a  lattice  structure  to  realize  linear  phase  transfer 
functions  and  develop  an  adaptive  algorithm  for  continuously  updating  the  lattice 
reflection  coefficients.  The  lattice  structure  is  considered  because  of  its  superior  finite 
wordlength  performance  compared  to  transversal  structures.  The  adaptive  lattice 
algorithm  developed  in  this  thesis  has  been  applied  to  estimate  the  sinusoidal 
frequencies  as  part  of  Prony's  method.  Results  of  computer  simulation  supporting  the 
theory  are  reported. 
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I.  INTRODUCTION 

Every  finite-impulse  response  (FIR)  filter  has  two  distinct  properties  LRef.  5j  , 
first,  it  is  always  stable;  second,  if  it  is  not  causal,  it  can  always  be  made  to  be  causal 
by  introducing  a  finite  delay.  FIR  filters  can  be  designed  so  that  their  frequency 
responses  have  an  exactly  linear  phase  characteristic.  FIR  filters  with  linear  phase  are 
important  in  applications  like  speech  processing  and  data  transmission,  because  in 
these  applications  a  nonlinear  phase  filter  is  harmful.  The  symmetry  property  of  the 
linear  phase  FIR  filters,  however,  helps  reduce  the  number  of  coefficients  by  nearly  one 
half  resulting  in  substantial  computational  savings. 

The  various  filter  realizations,  or  structures  that  are  frequently  considered  are  the 
direct  form,  the  cascade  form,  the  parallel  form,  and  the  lattice  form.  The  lattice  form 
realization  is  of  particular  interest  because  of  its  superior  numerical  performance,  and 
modularity  in  the  structure.  The  operation  of  a  lattice  filter  is  completely  described  by 
specifying  the  sequence  of  reflection  coefficients  that  characterize  the  individual  stages 
of  the  filter. 

Conventionally  an  adaptive  filter  is  composed  of  a  tapped  delay  line  or 
transversal  structure  with  adjustable  coefficients  or  weights  and  an  adaptive  algorithm 
which  updates  the  coefficients  continuously  based  on  some  performance  criterion.  The 
design  of  a  fixed  coefficient  filter  is  based  on  the  prior  knowledge  of  both  signal  and 
noise.  Adaptive  filters,  on  the  other  hand,  have  the  ability  to  adjust  their  own 
parameters  automatically,  and  their  design  requires  little,  or  no  a  priori  knowledge  of 
signal  or  noise  characteristics  [Ref.  14].  However,  the  designer  has  to  choose  the  order 
of  the  filter  and  the  type  of  the  algorithm.  Also  the  adaptive  filter  usually  requires  a 
large  initial  transient  time  (i.e.,  the  initial  filter  convergence  period). 

The  least-mean-square  fLMS)  adaptive  algorithm  minimizes  the  mean  square 
error  C!k.  oy  recursively  altering  the  filter  coeificient  vector  ?(k)  at  each  sampling 
mstant.  The  anginal  ’.Vihrov/- Hoff  LMS  algorithm  is  5(k-  i)  =  5(k)  — hjK\k)£(k/, 
where  n  is  a  convergence  factor  controlling  stability  and  rate  of  adaptation  [Ref.  15]. 
The  algorithm  is  based  on  the  method  of  steepest  descent,  moving  b(k)  in  proportion 
to  the  instantaneous  gradient  estimate  of  the  mean  square  error.  The  successive 
orthogonalization  provided  by  the  lattice  offers  advantages  in  adaptive  convergence 
rate  which  cannot  be  achieved  with  tapped  delay  lines  (Ref.  17,18]. 
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Recently,  Prony's  method  has  been  applied  to  the  estimation  of  spectral  lines  in 
noise  [Ref.  11].  It  has  been  shown  that  the  Prony's  method  yields  better  spectral 
estimates  than  a  companion  spectral  estimation  technique,  called  Pisarenko's  method 
(Ref.  10].  Also  it  has  been  established  that  the  filter  structure  involved  in  Prony's 
method  has  a  linear  phase  property  which  is  not  the  case  for  Pisarenko  s  method  [Ref 

HI- 

The  thesis  investigates  the  application  of  lattice  structures  in  Prony's  method  of 
spectral  line  estimation.  The  complete  solution  to  Prony's  method  consists  of  three 
steps:  (i)  representing  a  given  process  of  M  sinusoids  in  noise  in  terms  of  complex 
exponentials,  (ii)  finding  roots  of  a  symmetric  polynomial,  and  (iii)  estimating  the 
frequency,  phase  and  amplitude  information.  In  this  thesis,  however,  we  have 
emphasized  only  the  frequency  estimation  problem.  Here  we  derive  an  adaptive  lattice 
structure  to  realize  linear  phase  transfer  functions  which  will  be  used  to  estimate  the 
sinusoidal  frequencies  as  part  of  Prony's  method.  The  scope  of  the  thesis  consists  of 
obtaining  a  linear  phase  lattice  structure,  developing  a  least  mean  square  (LMS)  error 
based  adaptive  algorithm,  and  testing  the  lattice  structure  and  the  algorithm  by  means 
of  computer  simulation. 

The  thesis  is  organized  as  follows.  In  Chapter  II,  we  present  a  brief  review  of 
Prony's  method  of  representing  a  given  process  in  terms  of  a  set  of  complex 
exponentials,  and  then  address  the  problem  of  estimating  spectral  lines  using  this 
method.  We  show  that  the  original  Prony's  method  has  to  be  modified  slightly  in 
order  to  apply  it  to  the  spectral  line  estimation  problems.  In  Chapter  III,  we  discuss 
the  basic  concepts  of  the  linear  phase  FIR  filter  and  the  related  lattice  structure.  We 
show  three  examples  of  obtaining  the  lattice  structures  for  both  linear  phase  and  non¬ 
linear  phase  FIR  transfer  functions  (Appendix  A).  In  Chapter  IV,  we  briefly  discuss 
the  least-mean-square  (LMS)  adaptive  algorithm  that  results  from  attempting  to 
minimize  the  mean  square  error  and  present  a  summary'  of  some  LMS  algorithms  for 
lattice  that  have  been  reported  in  the  literature.  Also  included  are  some  computer 
simulation  results,  in  Chapter  V.  ve  present  'he  derivation  of  i  new  LMS  based  R 
viwUpti’rS  iiiiiiCw  tiiij cn tinci  extend  cC  me  iineur  pimsc  ense  es  wen.  new 
algorithm  is  then  used  in  the  estimation  of  spectral  lines  in  white  noise.  Results  of 
simulation  are  included. 
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II.  PRONY'S  SPECTRAL  LINE  ESTIMATION 

A.  INTRODUCTION 

In  its  original  form,  Prony's  method  analyzes  processes  involving  simple,  or 
damped  harmonics  using  complex  exponential  functions  as  coordinate  functions.  On 
the  other  hand,  the  well  known  Fourier  analysis  consists  of  representing  a  given 
process  in  terms  of  a  set  of  sine  and  cosine  functions.  Recently,  Prony's  method  has 
been  applied  to  the  estimation  of  spectral  lines  in  noise  [Ref.  11).  It  has  been  shown 
that  the  Prony's  method  yields  better  spectral  estimates  than  a  companion  spectral 
estimation  technique,  called  Pisarenko's  method,  in  terms  of  bias,  spurious  responses, 
and  the  computational  complexity  [Ref.  10]. 

In  this  chapter,  we  present  a  brief  review  of  the  Prony's  method  of  representing  a 
given  process  in  terms  of  a  set  of  complex  exponentials,  and  then  address  the  problem 
of  estimating  spectral  lines  using  this  method.  We  show  that  original  Prony's  method 
has  to  be  modified  slightly  in  order  to  apply  it  to  the  spectral  line  estimation  problems. 

B.  PRONY'S  METHOD 

Consider  that  the  given  process,  x(k),  consists  of  n  distinct  sinusoids,  then  x(k) 
can  be  approximated  by  an  expression  of  the  form 

x(k)  =.  i  [A.  coscoA  +  Bj  sintOjk]  (2.1) 

where  the  co/s  are  sinusoidal  frequencies.  The  above  approximation  can  be  considered 
as  a  special  case  of  an  exponential  approximation  given  by  [Ref.  29] 


x(k)  =  _!c  daiK  (2.2) 

where  m  =  2n  and  the  a/s  are  identified  as  C0j  and  -cOj.  The  values  of  (0;  can  be 
estimated,  by  Prony's  method,  assuming  that  the  data  are  known  at  k  =  0,  1,  .. ,  N-l. 
From  Eq.(2,2),  we  can  obtain  the  following  set  of  equalities 
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Cj  <Jal  +  c2  e^a2  +  ...  +  cm  eJam  -  x(l) 


C,  cJ2al  +  C,  ei2a2  +  ...  +  C  ei2an 
i  i  m 


Cj  eKN-Daj  +  q  eKN-Da2  +  ...  +  cm  eKN’1>am  ~  x(N-l) 

Now  wc  have  a  set  of  N  equations  with  2m  unknowns,  namely,  C  and  a.(  (i-  1,  2 . 

m)  which  can  be  solved  exactly  if  N  *  2m,  or  approximately  by  the  method  of  least 
squares  if  N  >  2m.  Also  note  that  the  N  equations  are  nonlinear  in  the  exponential 
terms  elai.  Let  ejai,  i*  1,  2,  ...,  m,  be  the  roots  of  the  equation  [Ref.  29]. 


eima  +  aiCj(m-l)a  +•  a2e^m'2^a  +  ...  +•  a^ei3  +■  am  -  0  (2.4) 

In  order  to  determine  the  coefficients  a.  (i  *  l,  2, ....  m),  we  multiply  the  first  equation 

in  Eq.{2.3)  by  am,  the  second  equation  by  am,j . the  mth  equation  by  a,  and  the 

(m+  l)1*1  equation  by  1,  and  add  the  results.  In  this  way,  we  can  obtain  the  (N-m) 
linear  equations.  From  Eq.(2.2),  we  can  obtain  the  following  set  of  equalities 

x(m)  +  x(m-l)a1  +  x(m-2)  a2  +  ...  +  x(0)am  *  0 
x(m-M)  +  xOm)^  +  x(m-l)a2  +  ...  +  x(I)am  *  0 


x(NT-l)  +  x(N-2)  al  +  x(N-3)  a2  +  ...  +  x(N'-m-l)  a  =  0 


Since  the  data  x(k),  k  =  0,  1,  2, N-l,  are  known,  this  set  of  equations  can  be  solved 
for  the  m  o's  if  N£2m.  After  the  a  s  are  determined,  eiai,  i  =  1,  2,  ....  m,  are  found 
from  Eq.(2.4).  The  set  of  equations  in  Eq.(2.3)  then  becomes  linear  in  terms  of  Cj,  i  = 
1,  2,  ....  m.  The  C's  can  be  determined  from  the  first  m  equations,  or  determined 
approximately  by  applying  the  least  squares  method  to  the  entire  set  of  equations.  The 
nonlinearity  associated  in  finding  elai,  which  is  related  to  the  frequencies  joy  and  -jto.,  is 
concentrated  in  Eq.(2.4).  The  above  procedure  is  known  as  Prony's  method. 

Now,  in  the  case  of  Eq.(2.1),  since  the  roots  of  Eq.(2.4)  occur  in  reciprocal  pairs, 
Eq.(2.4)  must  be  invariant  under  the  substitution  of  e'iai  for  eiai,  so  that  we  must  have 
a2n~  a(2n-l)”ai’  •"*  an+l“an-r  Thus  Eq.(2.4)  becomes 

ei2n<0  +  +  ...  +  an.j  +  an  eina> 

+  an^  ei(n*l)®)  +  ...  +  aj  J40  +  1  =0 

or 

ejnto  { (ejnco  +  g-jntO)  +  ^  (ej(n-l )co  +  g-jCn-l)^  +  _  + 

an-l  (ei<°  +  e'i<0)  +  «n  ]  =  0 

since  eJn<0  ¥"  0,  we  have 

2  cosnco  +  2aj  cos(n-l)o)  +  ...  +  2«n.j  cosco  +  an  =  0  (2.6) 

Now  noting  that  m=2n,  and  applying  the  above  symmetry  of  a  s  to  Eq.i2.5i  yields 
[Ref.  li] 


(x(0)  +  x(2n)}  +  {x(l)  +  x(2n-l)}  aj  +  ...  +  {x(n-l)  +  x(n+l)}  on.j 


+  x(n)  an  -  0 


{x(l)  +  x(2n+l)}  +  {x(2)  +  x(2n)}  aj  +  ...  +  {x(n)  +  x(n+2)}  an.j 


+  x(n+ 1)  an  -  0 


(2.7) 


{x(N-2n-l)  +  x(N-l)}  +  (x(N-2n)  +  x(N-2)}  a{  +  ...  +  fx(N-n-2) 
+  x(N-n)}  an_j  +  x(N-n-l)an  *  0 


Eq.(2.7)  consists  of  a  set  of  N-2n  equations  and  n  unknowns.  This  set  can  be  solved 
exactly  for  the  a's  if  N  *  3n,  or  solved  approximately  by  the  least  squares  method  if 
N  >  3n,  and  then  the  to's  are  determined  from  Eq.(2.6). 

C.  ESTIMATION  OF  FREQUENCY,  AMPLITUDE  AND  PHASE 

If  a  process  under  measurement  contains  an  unknown  number  of  sinusoids  of 
unknown  frequencies  and  amplitudes,  a  variant  of  Prony's  method  can  be  used  to 
determine  the  number  of  sinusoids  and  their  associated  frequencies  and  amplitudes.  As 
noted  above.  Prony's  method  is  a  technique  for  modeling  data  of  equally  spaced 
samples  by  a  linear  combination  of  exponentials.  An  exponential  curve  having  Yt 
exponentials  terms  can  be  determined  from  the  2M  data  measurements.  Each 
exponential  term  Aje3^  has  two  parameters  —  an  amplitude  Aj  and  an  exponent  a4 
where  a.  can  be  real  or  imaginary.  For  the  case  where  only  an  approximate  fit  with  M 
exponentials  to  a  data  set  ofN  samples  is  desired,  such  that  N>2M,  a  least  squares 
estimation  procedure  is  used. 


t  A. 


The  model  assumed  is  a  set  of  M  exponentials  of  arbitrary  amplitude,  phase, 
frequency  and  damping  factor.  A  process  consisting  of  M  real  undamped  (a  is 
imaginary)  sinusoids  can  be  expressed  as 


x(k)  -  j!  (a.z.k  +  ajVk)  (2.8) 

=  ^  A.  cos(2rrfkT  +  0.) 
i=l  ‘  1  1 

with  a.  *  (Aj/2)e^i  and  z.t  =  e^27t^,  where  A;  is  the  amplitude,  f.  is  the  frequency  and  0; 
is  the  phase  of  the  ic^  sinusoid,  respectively,  and  T  is  the  sampling  interval. 

Finding  {A.,  fj)  and  M  that  minimize  the  squared  error  is  a  difficult  nonlinear 
least  squares  problem.  Prony's  method  solves  two  sequential  sets  of  linear  equations 
with  an  intermediate  polynomial  rooting  step  that  concentrates  the  nonlinearity  of  the 
problem  in  the  polynomial  rooting  procedure  (similar  to  Eq.(2.4)). 

Define  the  polynomial,  B(z),  which  has  Z;  and  z*  as  its  roots,  given  by 


B(z)  -  n  (z  -  zXz  -  z .*)  -  ^  b.  z2M'i  =  0  (2.9) 

i=l  1  1  j  =  0  > 

with  bQ=  1  and  the  b-  being  real  coefficients.  Since  the  roots  are  of  unit  modulus  occur 
in  complex  conjugate  pairs,  then  Eq.(2.9)  must  be  invariant  under  the  substitution  z'* 
for  z.  Therefore,  Eq.(2.9)  can  be  written  as 


z2M  B(l/z) 


-  Z2M  b.  zi-2M 


i-o 


b.  z> 
=  0  J 


(2.10) 


Comparing  Eqs.(2.9)  and  (2.10),  we  may  conclude  that  bj  =  b2j^.j  for  j  =  0,1,  ...  ,M, 
with  bQ  =  b2^  =  1.  Thus  the  requirement  for  complex  conjugate  root  pairs  of  unit 
modulus  is  implemented  by  constraining  the  polynomial  coefficients  to  be  symmetric 
about  the  center  element.  Hence  the  realization  of  B(z)  is  a  linear  phase  FIR  filter. 
Based  on  order  2M,  a  linear  prediction  error  can  be  written  as 
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I 


$ 


(2.11) 


£(k)  -  f  b.  [x(k+j)  +  x(2M-k+j)] 
1-0  1 


which  reduces  the  number  of  coefficients  required  by  one-half. 

The  coefficients  b1(  b2,  ...  .  bm  are  determined  in  a  least  squares  fashion  by 
minimizing  the  total  squared  error. 


N-2M  o 
E  -  L  £“(k) 
i  =  0 


(2-12) 


which  yields  the  normal  equations 


e  b„ ,  ¥M- 


.  - „  \  [  X  [x(2M-k+i)  +  x(k+i)][x(2M-j  +  i)  +  x(j  +  i)]]  -  0 
]  =  0  K  1=0 


(2.13) 


for  k=  1 . M 

This  equation  can  be  solved  recursively  for  increasing  order  M. 

From  the  estimated  {b.}  values,  the  { z .}  are  determined  using  Eq.(2.9).  This  gives 


the  frequency  estimates. 


f.  =  tan1  [Im(z.)/Re(Zi)]  /  2ttT 


(2.14) 


To  obtain  the  {a.}  a  second  set  of  normal  equations  is  solved. 


N 


I  fa.  £  zj  z.l 
i-l  ;  ■  =  0  *  ! 


*  N  *:  *: 

t  Z,.  3  Z;  J  J 

•  .  =  0  1 


N 


X  0(2  Re  zkl)  x(j) 


(2.15) 


for  k-0,1, ...  ,M 
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III.  LINEAR  PHASE  FIR  FILTERING 


A.  INTRODUCTION 

Every  finite-impulse  response  (FIR)  filter  has  two  distinct  properties  [Ref.  5]  ; 
first,  it  is  always  stable;  second,  if  it  is  not  causal,  it  can  always  be  made  to  be  causal 
by  introducing  a  finite  delay. 

FIR  filters  can  be  designed  so  that  their  frequency  responses  have  an  exactly 
Linear  phase  characteristic.  FIR  filters  with  linear  phase  are  important  in  applications 
like  speech  processing  and  data  transmission,  because  in  these  applications  a  nonlinear 
phase  filter  is  harmful. 

The  impulse  response  sequence  of  a  linear  phase  FIR  filter  exhibits  a  kind  of 
symmetry,  e.g.,  h(n)  *  h(N-l-n)  for  n  =  0,  1,  ....  (N/2)-l  (assuming  that  N  is  even). 
In  general,  FIR  filters  require  a  large  number  of  coefficients  to  adequately  meet  with 
the  required  filter  specifications.  The  symmetry  property  of  the  linear  phase  FIR 
filters,  however,  helps  reduce  the  number  of  coefficients  by  nearly  one  half  resulting  in 
substantial  computational  savings. 

The  various  filter  realizations,  or  structures  that  are  frequently  considered  are  the 
direct  form,  the  cascade  form,  the  parallel  form,  and  the  lattice  form.  The  lattice  form 
realization  is  of  particular  interest  because  of  its  superior  numerical  performance,  and 
modularity  in  the  structure.  Lattice  realizations  have  been  successfully  utilized  in 
filtering  and  spectral  analysis,  and  in  modeling  of  some  physical  process  like  speech, 
geophysical  data  etc.  [Ref  8],  The  operation  of  a  lattice  filter  is  completely  described 
by  specifying  the  sequence  of  reflection  coefficients  that  characterize  the  individual 
stages  of  the  filter. 

In  this  chapter,  we  will  discuss  the  basic  concepts  of  the  linear  phase  FIR  filter 
and  the  related  lattice  structure.  And  finally,  we  will  show  three  examples  of  obtaining 

the  lattice  structures  realizing  for  both  phase  and  non-linear  phase  FIR  transfer 
functions  •  Appendix  A  . 

B.  LINEAR  PHASE  FIR  FILTERS  [REF.  3] 

Let  (h(n)}  be  a  causal  finite  duration  sequence  defined  over  the  interval 
O^n^N-1  having  the  linear  phase  symmetry  property  given  by  h(n)  =  h(N-l-n).  The 
z  transform  of  (h(n)}  can  be  written  as 


k  «  *  »  >  o. 


H(z)  *  £  h(n)z‘n 
n  =  0 


If  N  is  even,  then  Eq.(3.1)  can  be  rewritten  as 


<NP;'  h(n)z"n  + 


H(z)  “  £  h(n)z‘n  +  X 

n=0  n=N/2 


Now  applying  the  symmetry  property  of  h(n)  in  the  second  term  on  the  right  hand  side 
yields 


b(N/£1 


r»  +(NgH 


H(z)  -  X  h(n)z*n  +  X  h(N-l-n)z 
n=0  n=  0 

which  can  be  simplified  to 


H(z)  -  h(n)  {z-n  +  z-^'-1*0)} 

n  =  0 


If  N  is  odd,  then  Eq.(3.1)  becomes 


-(N-l-n) 


H(z)  «  I<N'^2H  h(n)[z'n  +  z'(N-'-">l  +  h(— )  z-t(N-l)/2] 
n=0  2 


If  we  evaluate  Eqs.(3.2)  and  (3.3)  for  z=  e^40  ,  we  obtain  the  discrete  Fourier  transform 
of  the  filter  sequence  h(n)  defined  as 


H(ej®>  -  X  h(n)  e'j*011 
n  =  0 


For  the  case  when  N  is  even,  we  then  have  the  discrete  Fourier  transform 


XjTXk  '  X  .  Xj  •*>  V  ' 


(3.5) 


H(e 


j«)-e-j  «{(N-D/2} 


(N/?)-l 

[  IT  2h(n)cos(  w{n-(N-l)/2)}]) 
n=*  0 


and  for  N  odd,  we  obtain 


H(ej<0)=,e'i<°^N’1)/2)[h{(N*1)'/2}  +  (NS3)/2  2h(n)  cos[co{n-(N-l )/2) ]].  (3.6) 

n  =  0 

In  both  cases  above,  the  sums  in  brackets  are  real,  implying  a  linear  phase  shift 
corresponding  to  a  delay  of  (N-l)/2  samples.  Figures  3.1  and  3.2  show  the  direct  form 
realizations  of  an  FIR  filter  with  linear  phase  for  both  N  even  and  N  odd  cases. 


Figure  3.1  Direct-form  realization  for  an  FIR  filter  with  linear  phase  (N  =  even). 


Figure  3.2  Direct-form  realization  for  an  FIR  filter  with  linear  phase  (N *  odd). 


C.  LATTICE  FILTERS  (REF.  4) 

The  basic  single  section  lattice  structure  is  shown  in  Figure  3.3  where  x(n)  is  the 
x(n)  and  fj(n)  and  gj(n)  are  generally  known  as  the  forward  and  backward  prediction 
errors,  respectively.  The  defining  equations  of  the  lattice  are  given  by 


f0(n)  *  go(n)  *  x(n) 
f,(n)  -  f0(n)  +  kjg0(n-l) 

g  1 1  n )  *  kjFji.n)  -  g^n-l) 


where  kt  is  called  the  lattice  reflection  coefficient.  If  another  section  is  cascaded  with 


the  first  one  the  resulting  equations  for  the  next  order  prediction  errors  are  given  by 


s 


i*.  s* * r  *" 


* *»  A  >  I'a  l'.  •**  »«.  |»,  I*.  »«  -u  ■-*  .1  t .1  ,#  >■*  k 


f2(n)  -  f^n)  +  kjg^n-l) 
g2(n)  -  k/^n)  +  gjCn-l) 


Substituting  for  fj(n)  and  gj(n-l)  from  Eq.(3.7)  yields 


f2(n)  *■  x(n)  +  (kj  +  ktk-.lxfn-l)  +  k,x(n-2) 


f2(n)  =  x(n)  +  b2  jx(n-l)  +  b2  2x(n_2) 


g2(n)“  k2x(n)+  {kj  +  kjk2}x(n-l)  +  x(n-2) 


g2(n)  *  b2i2x(n)  +  b21x(n-l)  +  x(n-2) 


where  b2  j  =  kj(  1  +  k2)  and  b2  2 -  k2- 


Thus,  by  induction,  for  a  cascade  connection  of  M  lattice  sections  we  have  the  general 


expressions 


fM<n>  “i  Jb  bM,ix(n-i} 


gM(°)  =4  bM,M-ix^n‘*> 


(3.10) 


where  q=  1.  Eq.(3.I0)  represents  the  FIR  filter  type  equations  to  obtain  the  Mth 
order  forward  and  backward  prediction  errors.  Now  taking  the  z-transform  of 
Eq.(3.10)  yields 


FM<Z>  "  ill,  bM,iz"  x« 
°M<Z>  bM,M.iz'‘  x<z) 


(3.11) 

(3.12) 


For  an  unit  impulse  input,  i.e.,  when  X(z)=  1,  F^(z)  and  Gyj(z)  represent  the  forward 
and  backward  prediction  error  filter  transfer  functions,  respectively,  and 


Gm<z>  -  z'Mfm<z'‘> 

l  _ 


(3.13) 


kM“  bM,M 


(3.14) 


In  fact,  the  FIR  like  prediction  error  filter  coefficients  can  be  iteratively  calculated 
starting  from  Eq.(3.14).  The  algorithm,  to  calculate  the  filter  coefficients  bu  also 

f  v  -V1  ** 

makes  use  of  the  well  known  lattice  property  that  an  order  lattice  contains  ail 

prediction  filters  of  orders  m^M.  Now  consider  the  mtb  section  (1  <m^  M)  in  the 
cascade  connection  of  M  lattice  sections  which  can  be  described  by  the  following 
equations: 


^m-l^z^  'r  kmz  ^m-l^2)  (3.15j 

Gm<z>  ‘  Wl®  +  z'‘Gm-l<z>  (3  >«) 

Eq.(3.16)  can  be  rewritten  as 


Gm.,(Z)  =  zGm(z)  -  zk^.^z)  (3.17) 

Substituting  Eq.(3.17)  into  Eq.(3. 15)  yields 


Fm<z>=  Fm-l(z>+ 


(3.18) 


t  Vi 

Therefore,  the  (m-l)ul  order  forward  prediction  error  transfer  function  can  be  written 
in  terms  of  the  mtb  order  forward  and  backward  prediction  error  transfer  functions  as 
follows: 


kmGm(z) 

k2 
*  m 

where  k^  =  '?T,  ;n  *  Recalling  Eq.i  3.15h  'hat  is. 
Gm(z)  =  z‘m  F^z'1) 


Fm-l(z)  = 


=  m 


Fm(z) 


1  - 


By  substituting  Eq.(3.20)  into  Eq.(3. 19)  we  find  that 


(3.19) 


(3.20) 


25 


»  •*  t.<  *  *>  i.1  «.•  >■>  i 


Fm-1^  " 


1  -  v- 


(3-21) 


Thus,  from  Eq.(3.21)  we  see  that  the  next  lower  order  polynomial  Fm.j(z)  can  be 
calculated  knowing  Fm(z).  Following  the  foregoing  procedure  we  can  find 
km_i  =  bm_i  from  Fm_j(z),  and  also  obtain  Fm„2(z).  This  way,  for  a  given 
order  polynomial  F^j(z),  we  can  determine  all  the  lattice  reflection  coefficients  1^, 
m=M,  M-l,  2,  1.  This  procedure  is  known  as  the  step-down  procedure  [Ref.  1], 
and  can  be  summarized  as  follows:  Let  the  given  order  polynomial  be 


FM<Z>  -  ;f0  'W1 


(3-22) 


and  by  replacing  M  with  m  we  have  an  m**1  order  the  polynomial  for  m  lattice  sections 


Fm<z>  -  ^0  W* 


(3.23) 


where  m*  M,  M-l, ....  2,  1.  We  now  define  another  polynomial  given  by 


Fm(z’')  -  ^0  Vf 


(3.24) 


As  we  step  through  the  procedure  from  m=  M  to  m=  1  Eq.(3.24)  can  be  expressed  as 


:  iZ'h  «  h  ,zm-i 

rrv  :*.f)  m.m-r 


th 

Define  km=  bm  m  and  obtain  the  m-l  order  polynomial  as  follows: 


SS 


»  *  »  •  •  »  *  *  *  *  •  *  ■ 
AW.V’i.S'N'V.S'A.V 


Wz> 


Fm(z)  *  kmZ'mFm  (Z'1  ) 


1  -  k2 


(3.26) 


Substituting  Eqs.(3.23)  and  (3.24)  into  Eq.(3.26)  yields 


jXvi.iz" 


j  J(;bnuz  1  ’  h™2  ' 


nbm-l,iz  1 


i  Jgbm,iz  1  *  *tmi^5m,m-iz  1 


(3.27) 


Then,  by  equating  the  coefficients  of  like  powers  of  z  on  both  sides  of  Eq.(3.27),  we 
obtain  the  computational  expression  for  the  (m-l)tb  order  polynomial  coefficients  as 


-’m- 1  ,i 


1  -  k2 


(3-28) 


where  m  =  M,  M-l,  ....  1  and  i  *  0,  1,  ....  m-1,  ^  “  bm  m  and  lkm'  <  1  *°r  a 
minimum  phase  polynomial  Fm(z). 

D.  LINEAR  PHASE  FIR  LATTICE  FILTER 

In  the  foregoing  we  considered  obtaining  a  lattice  structure  from  a  given  FIR 
transfer  function,  and  vice  versa.  In  what  follows  we  will  deal  with  a  special  case  of 
FIR  filtering,  namely,  the  linear  phase  FIR  filter.  Let  (h(n)}  be  a  causal  finite  duration 
linear  phase  sequence  defined  over  the  interval  O^n^N-1.  From  Eqs.(3.2)  and  (3.3). 
the  z  transform  of  him!  ran  he  vr.tten  as 


HN-1<z)  “  fm(z>  +  Z’D  gm(z> 


(3.29) 


S\S\VV'.V  '.VA,'.\VVA\S  .VA  .V.V.S  .V.V.V>.V.V1V 


^ kV.**  .'*•  k*» V*  .VV*\ 


Figure  3.4  Linear  Phase  FIR  Lattice  Filter. 


where  F^j(z)  ant*  G^j(z)  are  f°nvard  and  backward  filter  transfer  functions, 
respectively,  and  D  represents  the  number  of  unit  delays.  We  now  consider  the  N  odd 
and  N  even  cases  separately. 

iV  odd  case :  For  N  odd,  Eq.(3.29)  can  be  written  as 


HN-1<Z)  *  a0  +  a\z‘l  +  a2z'2  +  -  +  a(N-3)/2z‘(N'3)/2  +  a(N-l)/2z*(N*l) 


i') 


afN-3V2z 


*  ...  - 


aQ  +  ajz’1  +  a2z'2  +  ...  +  a{N.3)/2z'(N'3)/2  +  2  (1/2)  a(N.1)/2z*<N'1)/2 


+  a/x'  i\  ,z*<*+1)'2  +.  ...  +  a->z'^'3^  +  aiz*^'2^  +  anz*^'^ 


l(N-3),22 


28 


j 


Note  that  we  are  splitting  the  coefficient  l)/2  two  coefficients  of  value 

(l/2)a^_2)/2  eac^- 

HNt.|(z)  *  a0  +  ajZ*1  +  a2z*2  +  ...  +  a^.3^2z^N'3^2  0-30) 

+  0/2)  a(N-l)/2z*(N’1)/2  +  Z‘(N*1)/2  K]/2)a(N.1)/2 
+  a(N-3)/2z*1  +  -  +  a2z’(N°)/2  +  alZ-{N-3)'2 

+  a0z-(N-1>/2] 


HN-l(z)  "  Fm(z)  +  z*<N-D/2  Gm(z)  (3.31) 

where  the  order  of  the  polynomials  Fj^(z)  and  Gj^fz),  M  =  (N-l)/2,  and  the 
number  of  unit  delays,  D  =  (N-l)/2,  the  forward  transfer  function, 


FmW  =  a0  +  a,z-‘  +  ...  +  a(N.3)/2z-(N-3)/2  +  (l/2)a(N.])/2z-(N-')/2  (3.32) 

Fm(z'‘)  "  a0  +  a,z  +  ...  +  a(N.3)/2z(N'3)/2  +  (l/2)a(N.1)/2z(N-1)/2  (3.33) 


and  the  backward  transfer  function, 


►  v.t  (.l  l.t  l.l't.i  l.l  l  j'l.l  t.l  t  t.«~ Ll'l.i1 


iV  even  case:  For  N  even,  Eq.(3.29)  can  be  written  as 


Hn.i(z)  -  a0  +  ajZ*1  +  a2z-2  +  ...  +  a(N.2)/2z’<N-2)/2  +  a(N.2)/2z-N/2 
4-  ...  +■  a,z-(N-3>  4-  a.zW)  -  anz'^1) 


HN-l(z>  *  a0  +  alz'!  +  a2z*2  +  -  +  a(N-2)/2z' 


4*  z*N/'2  [a(N.2y2  4-  ...  4-  a 0z*{N*6)/2 


(N-2)/2 


(3.35) 


4-  a.z-(N-4)/2  4-  a0z-(N-2>/2  ] 


HNM(z)  =  Fm(z)  4-  z'N/2  Gm(z) 


(3.36) 


where  the  order  of  the  polynomials  F^(z)  and  G^j(z),  M  =  {(N/2)-l},  and  the 
number  of  unit  delays,  D  =  N/2,  the  forward  transfer  function. 


FM(z)  -  aQ  4-  ajz"1  4-  a2z‘2  4-  ...  4-  a^N.2^2z 


*(N-2)/2 


(3.37) 


FM(z  )  -  a0  +  ajz  +  a2z2  +  ...  +  a(N.2)/2z< 


(N-2)/2 


(3.38) 


and  the  backward  transfer  function. 


Gm(z)  -  z*<N"2)/2  Fm(z*1)  =  a(N.2)/2  +  ...  +  a2z' 

4-  ...  4-  aiz-(N‘4)/2  +a0z-(N-2)/2 


•(N-6)/2 


(3-39) 


Figure  3.4  shows  a  linear  phase  FIR  lattice  realization  where  the  filter  output  is 
obtained  by  adding  the  order  forward  prediction  error,  and  the  M**1  order  delayed 
(by  D  unit  delays)  backward  prediction  error.  Combining  the  discussion  in  this  section 
and  that  in  Section  (III.C),  we  can  summarize  the  algorithm  for  converting  a  given 
FIR  transfer  function  into  a  lattice  structure  as  follows: 

(i)  Find  if  N  is  even,  or  odd 

(ii)  For  N  odd:  M  -  (N-l)/2,  D  -  (N-l)/2 

(iii)  For  N  even:  M  *  {(N/2)-l),  D  »  N/2 

(iv)  From  F^z),  obtain  M  reflection  coefficients  as  discussed  in  Eqs.(3.22)  to  (3.26) 

(v)  Implement  the  lattice  as  shown  in  Figure  3.4 

E.  LATTICE  REALIZATION  OF  A  GENERAL  FIR  TRANSFER  FUNCTION 

In  the  foregoing  we  have  considered  a  polynomial  Fj^(z)  with  bp^  q“  1.  However,  in 
general  we  have  b^jQ^l.  In  this  section,  we  shall  modify  the  lattice  realization 
algorithm  presented  in  Section  (III.C)  to  suit  an  arbitrary  FIR  transfer  function  with 
0*1  22].  The  m1*1  order  polynomials  Fm(z)  and  Gm(z)  can  be  obtained 


from  £qs.(3.15)  and  (3.16): 

Fm<z)  "  smFm-l<z>  +  <3-40> 

^nfz)  “  kmFm-l(z)  +  sm  z  * <-'m- I  (z)  (3-41) 

where  the  coefficients  sm  **  bmQ  and  km  =  bm  m  are  recognized  as  the  reflection 
coefficients.  Eq.(3.41)  can  be  rewritten  as 

Gm-l(z)  "  s'm  z  Gm^  '  ^m  km  z  W2)  (3-42) 

Substituting  Eq.(3.42)  into  Eq.(3.40)  yields 

^m^z-  -  sm  Fm-1^  +  s  km^Gm^z^ '  kmFm-l^z^ 


r  H 

Therefore,  the  ^m-l/n  order  forward  prediction  error  transfer  function  can  oe 

wntten  in  terms  of  the  m^1  order  forward  and  backward  prediction  error  transfer 
functions  as  follows: 


From  Eq.(3.45),  we  obtain  the  computational  expression  for  the  (m-l)**1  order 
polynomial  coefficients  is 


h  .  =  I 

um-l,i  e2  .  ^2 


m 


m 


(3.46) 


where  m  -  M,  M-l . 1  and  i  -  0,  1 . m-1,  -  bmm,  sm  -  bm>0  and  sm>  ^ 

for  a  minimum  phase  polynomial  Fm(z).  It  may  be  noted  that  sm  =  1  for  m  *  1,  2, 
M-l.  This  indicates  that  we  have  only  one  s-coeflicient,  i.e.,  s^j,  which  requires  an 
extra  multiplication  in  the  M**1  lattice  section  as  shown  in  Figure  3.5. 
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IV.  LMS  ALGORITHM 


A.  INTRODUCTION 

Conventionally  an  adaptive  filter  is  composed  of  a  tapped  delay  line  or 
transversal  structure  with  adjustable  coefficients  or  weights  and  an  adaptive  algorithm 
which  updates  the  coefficients  continuously  based  on  some  performance  criterion. 

The  design  of  a  fixed  coefficient  filter  is  based  on  the  prior  knowledge  of  both 
signal  and  noise.  Adaptive  filters,  on  the  other  hand,  have  the  ability  to  adjust  their 
own  parameters  automatically,  and  their  design  requires  little,  or  no  a  priori  knowledge 
of  signal  or  noise  characteristics  [Ref.  14].  However,  the  designer  has  to  choose  the 
order  of  the  filter  and  the  type  of  the  algorithm.  Also  the  adaptive  filter  usually 
requires  a  large  initial  transient  time  (i.e.,  the  initial  filter  convergence  period). 

For  stationary  stochastic  inputs,  the  mean  square  error,  the  difference  between 
the  filter  output  and  an  externally  supplied  input  called  the  'desired  response",  is  a 
quadratic  function  of  the  filter  coefficients,  a  paraboloid  with  a  single  fixed  minimum 
point  that  can  be  sought  by  gradient  techniques  [Ref.  16j. 

In  the  previous  chapter  we  showed  that  the  operation  of  a  multistage  lattice  filter 
is  completely  described  by  specifying  the  sequence  of  reflection  coefficients  that 
characterize  the  individual  stages  of  the  filter.  In  this  chapter  we  briefly  discuss  the 
least-mean-square  (LMS)  adaptive  algorithm  that  results  from  attempting  to  minimize 
the  mean  square  error  and  present  a  summary  of  some  LMS  algorithms  for  lattice  that 
have  been  reported  in  the  literature.  Also  included  are  the  computer  simulation  results. 

The  least-mean-square  (LMS)  adaptive  algorithm  minimizes  the  mean  square 
error  e(k)  by  recursively  altering  the  filter  coefficient  vector  £(k)  at  each  sampling 
instant.  The  original  Widrow-Hoff  LMS  algorithm  is  B(k+  l)  =  5(k)  +  2ne(k)X(k), 
where  p  is  a  convergence  factor  controlling  stability  and  rate  of  adaptation  (Ref.  15). 
The  algorithm  is  based  on  the  method  of  steepest  descent,  moving  Vrt.-  in  •Proportion 
to  the  instantaneous  gradient  estimate  of  the  mean  square  error. 

When  the  filter  input  is  stationary,  the  backward  prediction  errors  are  orthogonal 
to  each  other,  with  the  result  that  the  successive  stages  of  the  lattice  filter  are 
decoupled  from  each  other  [Ref  8].  This  means  that  the  global  optimization  of  a 
multistage  lattice  filter  may  indeed  be  accomplished  as  a  sequence  of  local  optimization 


problems,  one  at  each  stage  of  the  lattice  filter.  Accordingly,  it  is  a  straightforward 
matter  to  increase  the  order  of  the  lattice  filter  by  simply  adding  one  or  more  stages 
without  affecting  the  earlier  design  computations.  The  successive  orthogonalization 
provided  by  the  lattice  offers  advantages  in  adaptive  convergence  rate  which  cannot  be 
achieved  with  tapped  delay  lines  (Ref.  17,18]. 

B.  .  SUMMARY  OF  THE  LMS  ALGORITHM 

The  LMS  algorithm  uses  an  estimate  of  the  gradient  of  the  mean  square  error  obtained 
from  the  adaptive  linear  combiner  which  is  a  combination  of  a  transversal  structure 
and  an  adder.  The  adaptive  linear  combiner  can  be  shown  in  two  basic  ways, 
depending  on  whether  the  input  is  available  in  parallel  form  (multiple  inputs),  or  in 
serial  form  (single  input)  [Ref.  6]. 


Figure  4.1  Parallel  Form. 


In  the  following  \vc  present  a  brief  derivation  of  the  LMS  algorithm.  Let  us  choose  the 
single  input  form,  then  the  filter  output  is  given  by 

y(k)  *  b0(k)x(k)  +  bj(k)x(k-l)  +  ...  +  bN.j(k)x(k-N+  1) 

N-l 

*  t  bj(k)  x(k-i) 

=  XT(k)  g(k) 

-  5  r( k )X[ k )  (4.n 

where  %  and  g  are  the  input  signal  vector  and  the  filter  coefficient  vector,  respectively 
and  (N-l)  is  the  order  of  the  filter.  The  input  signal  vector  X  and  the  filter  coefficient 
vector  g  arc  defined  as 
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Therefore,  the  output  y(k)  is  equal  to  the  inner  product  of  X(k)  and  B( k). 

The  error  e(k)  is  defined  as  the  difference  between  the  desired  response  d(k)  and 
the  actual  response  v(k), 


e(k)  =  d(k)  -  XT(k)B(k)  =  d(k)  -  ST(k)X(k)  (4.2) 

The  purpose  of  the  adaptive  algorithm  is  to  adjust  the  filter  coefficients  of  the  adaptive 
linear  combiner  to  minimize  the  mean-square  error.  A  general  expression  for  mean 
square  error  as  a  function  of  the  filter  coefficient  values,  assuming  that  the  input 
signals  and  the  desired  response  are  statistically  stationary  and  that  the  filter 
coefficients  are  fixed,  can  be  derived  in  the  following  manner  [Ref  6],  The  squared 
error  is 


e2(k)  =  d2(k)  -  2d(k)XT(k)B(k)  +  flT(k)X(k)XT(k)g(k)  (4.3) 

Taking  the  expected  value  of  both  sides  yields  the  mean  square  error, 


E[e2(k)j  =  E[d2(k)[  -  2  E[d(k)XT(k)]  B(k)  +  BT(k)  E[X(k)XT(k)J  B  (4.4) 

Defining  the  vector  P  as  the  cross-correlation  between  the  desired  response  and  the 
Input  vector,  we  have 

P  =  E[d(k)*(k)]  (4.5) 


Similarly  the  input  autocorrelation  matrix  R  is  defined  as 


(4.6) 


&  -  E[2(k)£T(k)] 

Thus  the  mean-square  error  can  be  expressed  as 

£(k)  *  E[e2(k)]  =  E[d2(k)]  -  2  ST  i(k)  +  §T(k)  fi  fi(k)  (4.7) 

The  gradient  V£(k)  of  the  mean  square  error  function  is  obtained  by  differentiating 
Eq.(4.7)  with  respect  to  the  filter  coefficient  vector  as  follows: 


V£(k)  = 


dE[e2(k)] 

5b0(k) 

<DE[e2(k)] 

5bv  i(k) 
k  ‘  1  . 


»  -2P  +2R  B( k) 


(4.S) 


The  LMS  algorithm  is  an  implementation  of  the  method  of  steepest  descent. 
According  to  this  method,  the  next  filter  coefficient  vector  is  equal  to  the  present  filter 
coefficient  vector  B(k)  plus  a  change  proportional  to  negative  of  the  gradient,  V£(k): 


fi(k+l)  -  g(k)-M£(k)  (4.9) 

The  parameter  H  is  the  factor  that  controls  stability  and  rate  of  convergence.  In  other 

words,  the  first  term  on  the  right  hand  side  consists  of  the  past  information  and  the 

second  term  represents  the  new,  or  updated  information. 

The  LMS  algorithm  estimates  an  instantaneous  gradient  in  a  crude  but  efficient 
«  ^ 

manner  ov  assuming  that  e"k).  the  square  of  a  singie  error  sample,  is  an  estimate  of 
the  mean-square  error  and  by  differentiating  e-,k)  with  respect  to  3(k).  flic  estimated 

gradient  is  given  by  the  following  expression 


A 


de^(k) 

<3b0(k) 


3e(k) 

5b0{k) 


Ve(k)  = 

=  2  e(k) 

• 

<3e^(k) 

5e(k) 

^bN-l^ 

The  estimated  gradient  components  are  related  to  the  partial  derivatives  of  the 
instantaneous  error  with  respect  to  the  filter  coefficient  components.  Thus  the 
expression  for  the  gradient  estimate  can  be  simplified  to 


Ve(k)  =  -2  e(k)  X(k)  (4.11) 

Using  this  estimate  in  place  of  the  true  gradient  in  Eq.(4.11)  yields  the  Widrow-HofT 
LMS  algorithm 

5(k+  1)  =  fl(k)  +  2  n  e(k)  X(k)  (4.12) 

Since  the  filter  coefficient  changes  at  each  iteration  are  based  on  imperfect  gradient 
estimates,  we  would  expect  the  adaptive  process  to  be  noisy,  that  is,  it  would  not 
follow  the  true  line  of  steepest  descent  on  the  performance  surface.  The  LMS  algorithm 
can  be  implemented  in  a  practical  system  without  squaring,  averaging,  or 
differentiation  and  is  elegant  in  its  simplicity  and  efficiency.  Each  component  of  the 
gradient  vector  is  obtained  from  a  single  data  sample  without  perturbing  the  filter 
coefficient  vector. 

C.  ADAPTIVE  LATTICE  ALGORITHMS 

The  parameters  :o  update  in  a  multistage  lattice  ire  its  reflection  confidents. 
Several  algorithms  have  been  proposed  in  me  literature  :br  upuuung  Tie  .uttice 

reflection  coefficients  [Ref.  7,8,17,18,20,23-28].  In  this  section  we  briefly  summarize 
some  of  those  algorithms.  Consider  a  lattice  filter  of  order  M.  For  stage  m  of  the 
filter,  the  flow  of  signals  is  described  by  the  following  pair  of  equations. 


Figure  4.3  A  Single  Stage  of  Lattice  Filter. 


fm<k>  ”  fm-l<k>  +  km,f  (4  ,3) 

Wk>  "  8m.l<k-!)  +  km,bfm-l<k) 

where  m=  1,2 . M 

For  stage  m,  the  forward  reflection  coefficient  is  denoted  by  km  j-  and  the  backward 
reflection  coefficient  is  denoted  by  km  ^ .  fm(k)  and  gm(k)  are  the  forward  prediction 
error  and  backward  prediction  error  of  stage  m  respectively.  Before  presenting  the 

adaptive  algorithms  for  the  lattice,  let  us  quickiy  summarize  some  non-adaptive 
retlection  coefficient  estimation  methods.  Later  on  we  can  obtain  the  adaptive 

updating  equations  as  approximations  to  these  methods. 


1 .  Non-adaptive  Methods 

When  the  lattice  coefficients  have  fixed,  non-adaptive  values,  several  methods 
have  been  proposed  for  computing  these  values  as  functions  of  the  correlation  statistics 
of  the  input.  Two  of  these  which  are  based  on  mean-square  error  (MSE)  minimization 
[Ref.2S]  are  given  beiow. 

Method  1  :  In  this  method,  we  choose  k^f  to  minimize  the  mean  forward 
prediction  error  power,  E[f2m(k)|  and  k^  ^  to  minimize  the  mean  backward  prediction 
error  power,  E[g2m(k)].  Here  we  give  the  final  equations  for  km  f  and  1^  ^  without 
going  into  the  details.  The  forward  and  backward  reflection  coefficients  are  given  by 


w 

Vb 


E[fm.|(k)gm.,(fc-l)l 

E[g2m.,(k.l» 

EtWk>gm.|(k-D] 


(4.14) 


where  both  km  f  and  km  ^  are  obtained  from  the  cross-correlations  between  the 
(m-l)^1  order  forward  and  backward  prediction  errors  normalized  by  respective 
prediction  error  powers. 


Method  2  :  For  a  single  channel  lattice  using  real  data  and  coefficients  we  can, 
however,  show  that  km  km  km-  Based  on  the  condition  that  km  £-=  km  ^  =  km, 
we  can  now  minimize  either  E[f2m(k)],  or  E[g2m(k)j  in  order  to  obtain  an  optimum  km. 
However,  it  seems  more  logical  to  minimize  the  sum  E[f2m(k)]  + E[g2m(k)]  as  suggested 
in  [Ref.  28].  The  resulting  reflection  coefficient  equation  in  its  final  form  is  given  by 


:  Eff^.;<;k)gm-!ik-[»i 


m- 


i  y  K- .  j  — 


or~  •  i  K.-  I 

°  :n-i 


where  we  have  normalized  the  cross-correlation  term  with  respect  to  the  sum  of  mean 
prediction  error  powers.  In  both  of  these  methods,  the  coefficients  at  stage  m  can  be 
computed  independently  of  those  following  that  stage,  Thus,  optimum  values  can  be 


computed  successively  along  the  structure  without  affecting  the  other  coefficients.  This 
phenomenon  lends  to  the  modular  nature  of  the  lattice  structure. 

2.  Adaptive  Methods 

An  instantaneous,  gradient-descent  adaptive  algorithm  minimizes  a  mean- 
square  error  criterion  can  be  derived  for  a  generalized  problem.  An  adaptive  lattice 
structure  is  suggested  by  incorporating  time  varying  coefficients  k^jfk)  and  km  ^(k)  in 
the  lattice  and  generating  algorithms  for  the  two  methods  described  above.  The 
resulting  procedures  are  : 

Method  l  :  Corresponding  to  method  1  of  non-adaptive  procedures  metioned 
above,  we  can  obtain  the  following  update  equations 


Wk+1> -  Wk> + 


ffm(k>  Sm-l^1)] 


(4.16) 


KmbOi+I)  ^mb*^  +  a2  av  ^m-l® 

a  m-l,g(k> 

where  p  is  an  adaptive  step  size  parameter,  X  is  a  positive  weighting  constant,  and  the 
forward  power  estimate  at  the  (m-l)^  stage,  f(k),  is 

«Wk>  sU2m-l/kl)  +  (1-X)  [g2m.l(k.l)] 
and  the  backward  power  estimate  at  the  (m-l)tk  stage,  <r2m.j  is 


Method  2  :  Eq.(4.15)  can  be  approximated  to  obtain  a  recursive  update  equation 
as  follows: 


k^k+i)  =  k^k)  +  p  [f^g^jfk-i)  +  rmA(Vgm(k)) 


(4.P) 


For  the  basic  structure  of  the  single-channel  lattice,  ie,  km  f=  km  k^,  the  reflection 
coefficients  or  filter  coefficients  km(k)  may  be  updated  using  a  modification  of  the 
LMS  algorithm  of  the  LMS  algorithm  [Ref.  17,20]. 

kn.(k+ 1)  =  l^flO  +  -i-  -  -  [fm(k)gm.,(k.l)  +  fm.,(k)gm(k)l  (4.18) 

a  m-lW 

where  is  the  power  estimate  at  the  (m-l)tk  stage.  Now  the  updated  power 

estimate  is 

<r2m.l(k)  =  X  ty^fk-l)  +  (l-M  te2m-l<k'1)  +  f2m-l<k)l  <4  l9> 

where  X  is  a  positive  weighting  constant  satisfying  the  criterion  0  <  X  5  1  then  controls 
the  bandwidth  of  the  filter  and  the  resulting  power  averaging  time.  A  power  estimate 
is  required  at  each  stage  in  the  lattice  due  to  the  fact  that  the  forward  and  reverse  error 
sequences  have  decreased  power  with  increased  stage  number. 

Method  3  :  A  third  successful  method  of  implementing  an  adaptive  algorithm  for 
FIR  lattice  structures  has  been  reported  by  Griffiths  [Ref.  17,18].  This  algorithm  has 
been  originally  discussed  for  a  noise  canceller  application.  The  form  of  the  lattice 
noise-cancelling  adaptive  filter  is  shown  in  Figure  4.4.  The  lattice  noise-cancelling 
adaptive  filter  consists  of  an  M  stage  linear  prediction  lattice  for  the  reference  signal 
xr(k)  together  with  a  set  of  tap  coefficients  Vm(k)  which  provide  the  noise-cancelling 
subtraction  paths.  Griffiths'  algorithm  is  briefly  presented  in  the  following.  The 
update  equations  for  km(k)  are  given  by 


Figure  4.4  Lattice  Form  Implementation  of  Noise-cancelling  Filter. 


Each  coefficient  Vn(k)  can  be  determined  independently  of  Vn(k)  for  n>m,  because  of 
orthogonalization  provided  by  the  lattice.  Thus  the  resulting  algorithm  is 


vm<k+  »  -  vm<k>  +  7/,,,  -  8m<k>l  <4-22> 

where  associated  power  measurement  y2m(k)  is 


y2m(k)  =  XY2m(k-l)  +  (l-X)[gm(k-l)gm(k-l)] 
f  h 

and  £m(k)  is  the  mul  stage  error  signal  as  shown  in  Figure  4.4. 


(4.23) 


D.  SIMULATION  RESULTS 


Figure  4.5  System  Identification  Experiment. 

For  the  purpose  of  computer  simulation,  let  us  consider  an  approximate 
algorithm  given  by 


km<k+  *>  *  km<k>  +  -  „  Aj-  [8m-l(k-')l  =(k)  («J) 

where 

<r2m(k)  =  X<r2m(k-1)  +  (1  -X)  g^k-l)  (*»-25) 
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The  configuration  used  in  the  simulation  is  a  system  identification  experiment  as  shown 
in  Figure  4.5.  The  fixed  filter  transfer  function  considered  is  given  by 
H(z)  -  1  -  0.89  z'1  +  0.25  z*2 

The  convergence  performance  of  the  LMS  algorithm  (as  given  by  Eq.(4.24))  can  be 
observed  by  plotting  the  error,  e(k).  versus  the  update  iteration,  k.  called  learning 
curves.  The  input  x(k)  to  both  fixed  and  adaptive  filters  is  a  white  noise  sequence  with 
zero  mean  and  unit  variance.  Figures  4.6  to  4.8  show  the  learning  curves  for  the  above 
example  where  we  have  used  three  different  values  for  the  adaptation  constant,  |i.  In 
the  next  chapter,  we  derive  a  new  algorithm  for  updating  the  reflection  coefficients, 
based  on  the  least  mean  square  principle  and  the  steepest  descent  algorithm. 
Improvement  in  convergence  speed  will  be  shown  using  the  new  algorithm. 


Figure  4.6  Learning  Curve  (ft  -  0.01). 


V.  ADAPTIVE  LINEAR  PHASE  LATTICE  ALGORITHM 


A.  INTRODUCTION 

In  Chapter  III  we  dealt  with  the  realization  of  fixed  coefficient  FIR  lattice  filter  with 
both  linear  and  nonlinear  phase  characteristics.  We  reviewed  the  basics  of  LMS 
algorithm  and  some  adaptive  lattice  realizations  reported  in  the  literature  in  the 
previous  chapter.  The  three  adaptive  lattice  algorithms  discussed  are  direct 
approximations  of  their  non-adaptive  counterparts.  None  of  them  estimates  a  gradient 
as  required  by  the  LMS  algorithm.  In  this  chapter  we  present  a  derivation  for  a  new 
LMS  based  FIR  adaptive  lattice  algorithm  and  extend  it  to  the  linear  phase  case  as 
well.  The  new  algorithm  is  then  used  in  the  estimation  of  spectral  lines  in  white  noise. 
Results  of  simulation  are  included. 


B. 


LMS  ALGORITHM  FOR  THE  FIR  LATTICE 

From  Eqs.(3.40)  and  (3.41)  the  FIR  lattice  equations  can  be  written  as 


Wk)  -  smfm-l(k>  +  iSngm-l^*1) 
Wk)  =  Wm-l^1)  +  kmWk) 


(5.1) 

(5.2) 


where  m=  1,2,.. .,M,  sm=  1  for  m^M  and  kj^  are  the  lattice  reflection  coefficients.  A 
realization  of  Eqs.(5.1)  and  (5.2)  is  shown  in  Figure  5.1.  The  lattice  input  is 
x(k)  =  gQ(k)  =  f()(k).  The  output  of  the  filter  is 


y(k)  =  fM(k) 

=  SM  ^M-l^  +  kM  §M-l^k'^ 


(5.3) 


Therefore,  :he  error,  e(k;,  is  given  by 


e(k)  =  d(k)  -  y(k) 

"  d(^  '  SM  fM-l^  *  kM  SM-l^’1) 


(5.4) 
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Figure  5.1  Adaptive  FIR  Lattice  Filter. 

where  d(k)  is  the  desired  signal.  The  objective  is  to  minimize  .he  mean  square  error 


J  -  E{c2(k)} 


(5.: 


The  gradient  of  J  with  respect  to  km  and  sm,  respectively,  are  given  by 


dJ  3\-(k) 

=  -2  Efc(k) 


c  k 


Ok: 


and 


*  •>.< 
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dJ 


ds 


M 


-2  E{e(k)  fM.j(k)} 


(5.7) 


Then  the  LMS  algorithm  can  be  formulated  as  follows: 


dJ 


kjdt-Kl)  «  kjCk)-^  — ) 


kj(k)  +  2  *ik  E{e(k)  (  -|Ll  )} 
kj(k)  +  2  pk  e(k)  ( -g-i  ) 


(5.8) 


and 


where  j=l,2,...,M.  The  next  step  is  to  estimate  the  gradient  vector  z(k).  For  this 
purpose,  define 


4>ij« 


15-12) 


and 


w  - 


Then  from  Eq.(5.3)  z(k)  is  given  by 


(5.13) 


z(k)  =  :(k) 

‘V1»J  (5  14) 

~  sM<k>  +  kM<k)  +  5M,j 

where  6^-  ==  (5k-yj/3kj).  Substituting  Eqs.(5.1)  and  (5.2)  to  Eqs.(5.12)  and  (5.13)  then 
we  have 


*y<k> 


Si(k) 


3k; 


ki(k> 


foi-lfr-t) 


3k: 


dk: 


fyfli)  -  Sj(k> 


SM <±R  .  wv&l® 


3k: 


ki(k> 


aki 


'.nererore. 


-  si<k)°i-l,j<k)  +  kiWi-ijtt-l)  +  SiAik-l)^- 

=  ^MjCk-l)  +  ki^MjW  + 

where 


(5.15) 

(5.16) 


Jh. 


(5.17) 


is  the  Kronecker  delta  function,  and  i  =  1,2,  ...  ,M  and  j  =  1,2,  ...  ,M.  If  i=0,  then 
Eqs.(5.15)  and  (5.16)  are 


®n  i00  - 


dx(k) 


(5.18) 


»0jW  -  ®o,j«  -  0 


(5.19) 


where  j  =  1 ,2,  ...  ,M.  Figure  5.2  shows  the  computation  of  gradient  elements  for  the 
adaptive  FIR  lattice  algorithm. 

From  the  Eqs.(5.15),  (5.16)  and  (5.19)  we  have 


oiti(k)  =  ^(k)  -  0 


(5.20) 


where  1  <i:S  j-1.  The  computations  of  gradient  elements  at  each  case  of  j  =  M,M-1, 
,1  are  as  follows: 

Case  j-  M  : 


•o.mW  =  *0lM<k)  =  0 

=  s;(k!(I)-.l.\pk»  -  k.( k)*P,_ j  vp. k- i )  -  gi.l(k-l)d-J<M 

Ti>M(k)  =  Si(k)'FM)M(k-l)  +  kifk^j.  1>M(k)  + 


where  i*  1,2,  ...  ,M.  From  Eq.(5.20),  every  gradient  element  equals  to  zero  except  final 
stage  elements  and  gradient  elements  of  the  final  stage  are 
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Figure  5.2  Computation  of  Gradient  Elements  for  the  Adaptive  FIR  Lattice  Algorithm. 


*M,M^  "  sM^k),^M-l,Ni(k*1)  +  kM<k)<I>M-l,M(k)  +  fM-l(k)6M,M 
Applying  the  Eqs.(5.17)  and  (5.20)  yields 
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1 


r, 


* 


r"< 


>3 


& 


3 


^M.mW  =  fM-lOO 


(5.21) 

(5.22) 


Case  j—  M-l  : 


'•'i.M-lW  -  ¥k)'*'i.l,M-l(k-1)  +  ki<k)®H,.M-l(k)  +  ri-l<k)«i,M-l 


where  i=  1,2,  ...  ,M.  Applying  the  Eq.(5.20),  last  two  stage  elements  are  considerable, 
if  i=  M-l,  then  we  have 


°M-l,M-lW  “  sM-l(k^M-2,M-l^k)  +  kM-l(^'PM-2,M-l(k‘1)  +  SM^'^M-l.M-l 

^M-l.M-l^  "  sM-lW'FM-2,M-l^k*1)  +  kM-l^<I>M-2,M-lW  +  fM-2(k)6M-l,M-l 
Applying  the  Eqs.(5.17)  and  (5.20)  yields 


°M-1,M-Uk)  “  gM-2^'1)  (5<23) 

TM-1,M-1^  =  fM-2(k)  (5-24) 

And  the  last  stage  terms  are 

=  ~  kMVk/FM-l,M-l(k*^  ~ 

TM,M-l(k)  =  sM(k)<FM-l,M-l(k'1)  +  kM<k^M-l,M-lW  +  fM-l<k)6M,M-l 


Using  the  Eqs.(5.17),  (5.23)  and  (5.24)  yields 


*  sM(k^M-2<k*1)  +  kM<k)  fM-2(k'^ 
^M.M-l^)  *  sM(k)fM-2^k'!)  +  kM(k)  8M-2(k_1) 

Finally, 

Case  j**  l  : 

®0,1«  -  *0,1®  -  0 


*U®  -  «i®®i-U®  +  ^*1-1,1®*)  +  *i-l®‘>*i.l 
*U®  "  *1®*!-!.!®')  +  ki<k)®i.1,l®  +  <i-l®*U 

where  i*  1,2, ...  ,M. 

The  LMS  algorithm  derived  in  the  foregoing  can  be  summarized  in  Table  1. 

Simulation  Results  : 

The  performance  of  the  algorithm  summarized  in  Table  1  has  been  observed  by 
computer  simulation.  The  configuration  used  is  the  system  identification  experiment  as 
discussed  in  Section  (IV.  D)  using  the  same  fixed  coefficient  filter.  The  learning  curves 
obtained  using  the  new  algorithm  are  shown  in  Figures  5.3  to  5.5.  Comparing  these 
with  the  learning  curves  in  Figures  4.6  to  4.8,  which  are  obtained  using  approximate 
algorithms,  we  observe  significantly  faster  convergence  rate  for  the  new  algorithm. 


(5.27) 

(5.28) 


(5.25) 

(5.26) 


»1<. 


t.iUH 


TABLE l 

LMS  ALGORITHM  FOR  THE  FIR  LATTICE 
Initialization: 

K(  0)  =o 

SM{0)  =  1 
Lattice : 

x(k)  =  fQ(k)  =  g0(k) 

=  Vk)fm-l<k>  + 

*m(k>  =  sm<k)gm-l(k‘1)  +  Vk>fm-l(k) 
m  =  1,  2,  .  .  .  ,  M 

y(k)  =  fM(k) 

Update  Equations: 

kj(k+l)  =  kj(k)  +  2  (*ik/<r2k.(k)]  e(k)  2j(k) 

SjkjI  k+  i )  =  sM(k)  +  2  [Ms/<T2s(k)]  e(k)  fM-1(k) 
where  j»k  and  fig  are  the  adaptation  constants, 

<x2k.(k)  =  X  <r2kj(k-l)  +  ( 1-X)  02M  j(k) 
and 

<r2s(k)  =  X  ff2s(k-l)  +  ( 1-X)  f2M.!(k) 

are  estimations  of  power  in  Zj(k)  and  fM-l(k), 
respectively  and  X  is  a  positive  weighting  constant, 
0<Xil. 

Gradient  Vector  Elements : 

*0,j(k)  =  ^0,j(k)  =  0 

<I>i  j(k)  =  si(k)4>i.1<j(k)  +  ki(k)'Fi.1/j(k-l) 

>  g,  ■. i ( k-1 )ot  , 

=  3i(k)'Fi.1>  ^k-l)  +  ki(k)<Di.1/ j(k) 


zH(k) 


■ 


’-f-T. 


2.5 
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Figure  5.3  Learning  Curve  (fi  =  0.1). 


Figure  5.5  Learning  Curve  (p  *  0.5). 


C.  LINEAR  PHASE  FIR  LATTICE  ALGORITHM 

Wc  now  extend  the  LMS  algorithm  derived  in  the  previous  section  to  the  linear 
phase  FIR  lattice  Filter. 

From  Eq.(3.29),  output  of  the  linear  phase  FIR  lattice  can  be  written  as 


y(k)  *  fM(k)  +  gM(k-D)  (5.29) 

The  lattice  input  is  x(k)  =  f()(k)=  gQ(k).  Substituting  Eqs.(5.1)  and  (5.2)  in  the  output 
equation  of  the  filter  yields 


y^k)  -  sM  fM.!(k)  +  kM  gM_ i<k- 1  >  +  sM  gjvi-l(k-D-l)  +  kM  fM-l*k‘D|5  30) 
“  SM  ifM-l(k)  +  8M-l(k-D-l)]  +  kM  ffM-l(k'D)  + 

and  the  error,  e(k),  is  given  by 
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Figure  5.6  Adaptive  Linear  Phase  FIR  Lattice  Filter. 
e(k)  -  d(k)  -  y(k) 

*  d(k)  -  sM  lfM.!(lt)  +  8M-l(k-D-l)l  *  kM  lfM*l(k’D)  +  SM-l(k'n<l 

where  d(k)  is  the  desired  signal.  A  schematic  of  the  linear  phase  FIR  lattice  realization 
is  shown  in  Figure  5.6.  From  Eq.(5.31)  we  see  that  y(k)  is  a  function  of  both  s^  and 
k^j.  But  f;vi_i(k)  and  g^j.^k-l)  are  functions  of  and  so  on.  Therefore,  in  order 
to  minimize  the  mean  square  error,  we  define  a  cost  function 


i  =  E  {  e2  (k)  ’ 


(5.32) 


and  a  gradient  element 


k«  fa  jAAti  I'ljt  I* 


fo(k) 

5ki 


(5.33) 


Then  following  the  treatment  in  Eqs.(5.6)-(5.ll)  the  LMS  algorithm  for  the  linear 
phase  lattice  can  be  written  as 


kj(k+  1)  -  k^k)  +  2  Hk  e(k)  z(k) 

sM(k^  1)  -  *M(k)  +  2  ms  e(k)  (fM.i(k)  +  gM.1(k-D-l)] 


(5.34) 


where  j*=  1,2,. ..,M.  From  Eqs.(5.29)  and  (5.30),  the  gradient  vector  z(k)  is  given  by 


(535) 

‘  SM»  +  +  kM<k>  I®M.l,j<k'D) 

+  TM.,,j(k-l)l  +  (fM.;(k-D)  -  gM.,(k.01«M  j 

where  ^(dk^/dkj),  and  <Pmj  and  *Fmj  are  obtained  on  the  same  lines  as  in 
Eqs.(5.15)-(5.17).  The  resulting  LMS  algorithm  is  summarized  in  Table  2. 

Simulation  Results  : 

The  performance  of  the  algorithm  summarized  in  Table  2  has  been  observed  by 
computer  simulation.  The  configuration  used  is  the  system  identification  experiment  as 
discussed  in  Section  (IV.  D).  We  have  used  two  different  fixed  coefficient  linear  phase 
FIR  filter  examples,  one  with  N  even  and  the  other  with  N  odd.  They  are 
H3(z)  -  0.154  +  0.462Z*1  +  0.462z'2  +  0.154z*3 
and 

H4(z)  =  0.15 -0.45Z'1  +  0.36z'2  -  0.45z'3  +  0.l5z‘4 

The  '.earning  curves  obtained  using  the  new  algorithm  are  shown  in  Figures  5.-  to  5.10. 
Figures  5."  and  5.3  are  obtained  learning  curves  with  linear  phase  FIR  transfer 

function  H3(z).  Figures  5.9  and  5.10  are  obtained  learning  curves  with  linear  phase 
FIR  transfer  function  H4(z). 


\y. 

‘.v. 

w, 

■  Sf 


TABLE  2 

LINEAR  PHASE  FIR  LATTICE  ALGORITHM 
Initialization: 

5(0)  =o 

*M<0>  =  1 

Lattice: 

x(k)  =  f0(k)  =  g0(  k) 

fm<k>  =  Vk>fm-l<k>  + 

9m<k>  =  +  Vk>fm-l<k> 

m  -  1,  2,  . . .  ,  M 

y(k)  =  sM  [fM_x(k)  +  gM-1(  k-D-1)  ]  +  kM  [  fM-1(  k-D) 
Update  Equations : 

kj(  k+1)  =  kj(k)  +  2  [nk/<r2k.(k)]  e(  k)  Zj(k) 
sM(k+l)  =  sM(  k)  +  2  (ns/(T2a(k)]  e(k)  (f^k) 

+  g,*.!*  k-D-1)] 

where  nk  and  n3  are  the  adaptation  constants, 
«r2kj(k)  =  X  <r2kj(k-l)  +  ( 1-X)  [0M  :j(k)  ^(k-D)] 

and 

<r2g(k)  =  X  <J2S(  k-1)  +  ( 1-X)  [fM-1(k)  +  g^k-D-l)] 

are  estimations  of  power  in  Zj(k)  and 
[  fM_l(  k)+gM_l(  k-D-l) ], respectively  and  X  is 
a  positive  weighting  constant,  0<X£1. 

Gradient  Vector  Elements: 

*o,j(k>  =  *o.j<k>  =  0 

<Pi(i(k)  =  si(k)<I>i.1  ,-(k)  *  ki(k)'Pi-1<  ,(k-D 

*  '^i-l(  k”1  )^i,  j 

Ti,j<k>  =  ^^i-l  #j(k-l)  +  ^(kjOi.i^k) 

+  fi-l<k)5i,j 
=  <J>M  ,<k)  +  ,(k-D) 


zH(k) 


Figure  5.8  Learning  Curve  (n  =  0.3). 
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Figure  5.9  Learning  Curve  (ji  =  0.03). 


ITERATIONS 


Figure  5.10  Learning  Curve  (fi  =  0.1). 


D.  SPECTRAL  ESTIMATION 

In  this  section,  we  extend  the  linear  phase  FIR  adaptive  lattice  algorithm  derived 
in  the  previous  section  to  the  estimation  of  spectral  lines  in  white  noise.  The  spectral 
estimation  problem  is  some  what  different  from  the  system  identification  experiment 
that  we  have  been  considering  so  far.  In  a  typical  spectral  estimation  problem,  we  are 
given  a  time  series  and  we  need  to  estimate  its  spectrum.  A  suitable  configuration  that 
is  frequently  used  for  this  purpose  is  shown  in  Figure  5.11.  Comparing  this  with 
Figure  4.5,  we*  notice  that  we  have  the  given  time  series  at  the  input  of  the  adaptive 
filter  and  the  filter  updates  its  coefficients  to  minimize  the  mean  square  output.  In 
doing  this,  we  assume  that  the  input  process  x(k)  has  been  generated  by  passing  a 
white  noise  sequence  through  a  filter,  say  I(z),  and  the  adaptive  filter  attempts  to 
realize  an  inverse  filter,  say  H(z)=l/I(z).  Considering  that  the  adaptive  filter  has 
sufficient  degrees  of  freedom,  the  output  of  the  filter  will  be  a  white  noise  sequence. 

The  adaptive  algorithm  summarized  in  Table  2  can  be  used  in  spectral  estimation 
after  appropriate  modification.  In  the  present  case  we  minimize  the  cost  function, 

J  =  E{y2(k)} 

rather  than  J  =  E{e2(k)}.  The  resulting  update  equations  can  be  shown  to  be 

kj(k+  1)  =  kj(k)  -  2  [Hk/<r\.(k)]  e(k)  Zj(k)  (5.36) 

sM(k+  1)  =  sM(k)  -  2  [*is/<r2s(k)l  e(k)  [^(k)  +  gM.1(k-D-l)]  (5.37) 

where  Zj(k),  <*k.(k)  and  <rs(k)  are  as  defined  in  Table  2. 

Using  the  adaptive  algorithm  in  Eqs.(5.36)  and  (5.37)  we  can  obtain  values  of  the 
reflection  coefficients  and  the  s^j  coefficient.  We  then  obtain  the  equivalent 

polynomial.  Bfz),  by  going  through  the  standard  step  up  procedure  given  by  [Ref.  !]: 


where  i~  l,2,...,m-l  and  m=  Finally  the  required  linear  phase  transfer 

function  is  obtained  as  follows: 


HN-l(z>  =  FM<z>  +  Z*D  gm(z> 

where  Fjy((z)  =  B(z)  and  G^jCz)®  z‘^B(z**).  The  spectrum  is  computed  as  follows: 


1 

S(0 - 

i  ho  +  e-i1(0  +  ...  +  hN.j  e-K11*1^  I2 

where  <o  =  2ttf,  and  f  is  the  normalized  frequency  (with  respect  to  the  sampling 
frequency)  in  the  range,  0  £  f<  0.5. 

Simulation  Results  : 

The  input  process  x(k)  consists  of  a  signal  in  noise,  given  by 
x(k)  *  s(k)  +  w(k) 

where  s(k)  may  be  a  single,  or  multiple  sinusoids  and  w(k)  is  a  zero  mean  unit  variance 
white  noise  sequence.  In  the  following  we  consider  several  examples  using  parameters 
ranging  from  a  single  sinusoid  to  4  sinusoids,  SNRs  from  30dB  to  lOdB,  and  filter 
order,  M,  from  2  to  30. 

Example  I:  We  consider  a  single  sinusoid  given  by 


x(k)  =  cos(2ttfk)  +  w(k)  (5.38) 

where  f=0.l5,  SNR=30dB. 

Figures  5.12.  5.13,  5.14,  and  5.15  are  plots  of  frequency  spectrum  with  different 
order  of  lattice  M  and  adaptation  constant  p. 

Figure  5.12  is  the  plot  of  the  case  M  =  2  and  n  =  0.015.  We  see  that  the  peak  is 
at  f=0.15  and  no  spurious  responses. 

Figure  5.13  shows  the  frequency  spectrum  of  M  =  4  and  p  =  0.01.  We  observe 
that  one  peak  is  at  f  =  0.15  and  a  spurious  response  whose  magnitude  is  about  -53dB  at 
f=0.37. 


Figure  5.14  shows  the  frequency  spectrum  of  M  *  10  and  fi  =  0.03.  There  arc  one 
peak  at  f=0.15  and  four  spurious  responses  at  f=0.05,  0.25,  0.35  and  0.45.  The  largest 
spurious  response  is  about  -15  dB  at  f=0.45. 

Figure  5.15  shows  the  frequency  spectrum  of  M  =  20  and  =  0.03.  There  arc  one 
peak  at  f=*0.15  and  nine  spurious  responses.  The  largest  spurious  response  is  about 
-27  dB  at  f=  0.125. 

Through  the  Figures  5.12,  5.13,  5.14,  and  5.15  we  can  determine  that  the  best 
model  order  to  detect  the  one  sinusoidal  signal  is  VI  =  2. 

Next,  we  fix  the  order  of  lattice  M.  adaptation  constant  u  and  iteration  number 
k.  Figures  5.16,  5.17  and  5.18  are  the  plots  of  frequency  spectrum  with  changing  the 
signal  to  noise  ratio  (SNR). 

Figure  5.16  is  the  plot  of  M  =  20,  n  =  0.05,  k=3000  and  SNR=30dB.  There  is  a 
peak  at  f*0.15  and  the  largest  spurious  response  is  about  -27dB  at  f=  0.125. 
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Figure  5.17  is  the  plot  of  M  =  20,  ft  =  0.05,  k=  oOOO  and  SNR  =  20dB.  The  largest 
spurious  response  at  f=  0.425  is  larger  than  the  desired  response.  At  this  case,  we 
cannot  estimate  the  frequency  of  sinusoid. 

Figure  5.18  is  the  plot  of  M  =  20,  p.  =  0.05,  k=3000  and  SNR=  lOdB.  Three  of 

the  spurious  responses  are  larger  then  the  desired  response.  The  desired  response  is 

about  -13dB. 

From  Figures  5.16,  5.17  and  5.18,  if  the  SNR  is  getting  smaller  then  estimating 
sinusoidal  frequency  is  more  difficult. 

Example  2:  Next,  to  estimate  two  sinusoidal  frequencies  in  noise,  set  the  input  x(k)  as 


x(k)  =  V?  {  cos(27Tf,k)  +  cos^Ttfjk) }  4-  w(k) 


(5.39) 


where  the  normalized  frequencies  of  signals  fj  =  0.15  and  f2=0.25  and  set  SNR=30dB. 

Figure  5.19  shows  the  frequency  spectrum  of  M  =  4  and  u  =  0.02.  There  are  two 
peaks  at  f ^  =  0.15  and  f2  =  0.25  and  no  spurious  responses. 

Figure  5.20  shows  the  frequency  spectrum  of  M  =  20  and  ji  =  0.022.  There  are 
two  peaks  at  fj  =  0.l5  and  f2  =  0.25,  and  eight  spurious  responses.  The  largest  spurious 
response  is  about  -20dB  at  f=  0.375. 

From  Figures  5.19  and  5.20,  the  best  model  order  to  estimate  the  frequencies  of 
two  sinusoidal  signals  is  M  =  4. 

Example  3\  Next,  to  estimate  four  sinusoidal  frequencies  in  noise,  set  the  input  x(k)  as 

x(k)  =  {  cos(2rtfjk)  +  cos(2rtf2k)  +  cos(27tf^k)  +  cos(27tf^k)  }  +  w(k)  (5.40) 

where  the  normalized  frequencies  of  signals  f ,  =  :\  =  d  ;5.  f-  =  > >. 25.  N  = ’\35.  inti 

set  SNR  =  30dB. 

Figure  5.21  shows  the  frequency  spectrum  of  M  =  8  and  fl  =  0.015.  There  are  four 
peaks  at  fj  =  0.07,  f2  =  0. 1 5,  ^  =  0.27  and  =  0.375,  and  no  spurious  responses. 

Figure  5.22  shows  the  frequency  spectrum  of  M  =  30  and  p.  =  0.0 14.  There  are 
four  peaks  at  fj  =  0.05,  f2  =  0.15,  ^  =  0.25  and  fj=0.35,  and  eleven  spurious  responses. 

The  largest  spurious  response  is  about  -23dB  at  f=0.45. 
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From  Figures  5.21  and  5.22,  the  best  model  order  to  estimate  the  frequencies  of 
four  sinusoidal  signals  is  M  =  8. 
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Figure  5.15  Frequency  Spectrum  (J  sinusoid.  M  =  20,  ji  =  0.03). 


Figure  5  16  Frequency  Spectrum  ( l  sinusoid,  M  *  20,  SNR=30dB). 
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Figure  5.17  Frequency  Spectrum  (t  sinusoid,  M  =  20,  SNR  =  20dB). 


Figure  5.18  Frequency  Spectrum  (1  sinusoid,  M  =  20,  SNR=  lOdB). 


Figure  5.19  Frequency  Spectrum  (2  sinusoids,  M  =  4,  [1  =  0.02). 


Figure  5.20  Frequency  Spectrum  (2  sinusoids,  M  =  20,  p  =  0.022). 
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Figure  5.21  Frequency  Spectrum  (4  sinusoids,  M  =  8,  fl  =  0.015). 


Figure  5.22  Frequency  Spectrum  (4  sinusoids,  M  *  30,  )i  =  0.014). 
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E.  SUMMARY,  CONCLUSION  AND  SUGGESTED  FUTURE  WORK 

The  thesis  investigates  the  application  of  lattice  structures  in  Prony's  method  of 
spectral  line  estimation.  The  complete  solution  to  Prony's  method  consists  of  three 
steps:  (i)  representing  a  given  process  of  M  sinusoids  in  noise  in  terms  of  complex 
exponentials,  (ii)  finding  roots  of  a  symmetric  polynomial,  and  (iii)  estimating  the 
frequency,  phase  and  amplitude  information.  In  this  thesis,  however,  we  have 
emphasized  only  the  frequency  estimation  problem. 

The  underlying  theory  of  Prony's  method  has  been  briefly  reviewed  in  Chapter  II. 
By  using  a  modified  Prony's  approach,  we  observed  that  the  frequency  estimation  (as 
part  of  Prony's  method)  requires  a  polynomial  with  a  linear  phase  property.  The 
basics  of  the  linear  phase  FIR  filter  and  the  lattice  structure  have  been  included  in 
Chapter  III.  Also  we  have  shown  that  a  linear  phase  FIR  transfer  function  can  be 
realized  from  an  FIR  lattice  using  D  number  of  additional  unit  delays  and  an  adder, 
where  D  is  determined  by  the  order  of  the  linear  phase  polynomial. 

The  principles  of  adaptive  filtering  and  the  LMS  algorithm  have  been  briefly 
studied  in  Chapter  IV.  We  have  derived  a  new  LMS  algorithm  for  the  FIR  non-linear 
phase  ans  linear  phase  transfer  functions  in  Chapter  V.  The  problem  spectral 
estimation  using  the  new  adaptive  algorithm  has  been  addressed,  and  the  simulation 
results  are  included. 

The  significant  outcome  of  the  work  is  the  derivation  of  an  LMS  type  algorithm 
for  a  linear  phase  lattice  structure  and  its  application  to  spectral  line  estimation  as  part 
of  Prony's  method.  As  has  been  shown,  the  new  algorithm  yielded  faster  convergence 
rate  compared  to  previously  reported  approximate  algorithms.  The  application  of  this 
algorithm  to  spectral  estimation  produced  high  spectral  resolution  as  illustrated  by  the 
simulation  results. 

However,  we  have  observed  spurious  spectral  responses,  when  the  model  order  is 
higher  that  the  required.  Also  we  have  observed  that  the  SNR  performance  of  the 
algorithm  needs  to  be  improved.  For  input  SNRs  less  than  about  10  dB.  the  algorithm 
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APPENDIX  A 

EXAMPLES  OF  OBTAINING  A  LATTICE  STRUCTURE  FOR  THE  FIR 

TRANSFER  FUNCTION 

Example  A.  I 

Obtaining  a  lattice  structure  for  the  general  FIR  transfer  function.  The  unit 
sample  response  of  a  FIR  transfer  function  is  given  by 
H3(z)  -  0.5  +  0.25Z*1  +  0.125z'2  +  0.0625z'3 

Solution :  Given  unit  sample  response  is 

H3(z)  «  0.5  +  0.25z'1  +■  0.1 25z'2  +  0.0625z*3  (A.l) 

Using  Eq.(3.22),  we  have  polynomial  for  3  sections 

H3^  "  b3,0  +  b3,lz'1  +  b3,2z'2  "  b3,32'3  (A-2) 

Comparing  Eqs.(A.l)  and  (A.2),  we  have 
b3,0  *  0  5 

b3.1  ”  0.25 

(A. 3) 

b32  =  0.125 


b3  3  -  0.0625 

Starting  with  m~  3  we  have  from  Eq.(3.46) 


*3  =  °3,3  =  °'0b25 
s3  *  b3,0  *  0  5 


(A. 4) 


Now,  we  need  to  generate  the  coefficients  for  H2(z)  and  from  Eq.(3.46) 
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-smbau  ~  *mbm,m-i 


Vu  - 

sm  ' 


And  for  m*  3  and  i-  0  we  have 


52,0 


Jib,1)'  k3b3  3 

«  2  _  v-2 


(A. 5) 


(0.5)(0.5)  -  (0.0625)(0.0625) 
(0.5)2  -  (0.0625)2 


(A.  6) 


Next,  for  m  =  3  and  i  =  1  we  get 


s3b3.1  '  k3b3.2 

— ri — 7~r 

s3  •  k3 


(0.5){0.25)  -  (0.0625)(0.125) 
((X5)2  (0.0625)2 


=  0.47619 


and  finally 


j3b3.,2  -  k3b3,l 

-  k3^ 


(A. 7) 


(0.5)(0.125)  -  (0.0625X0.75) 
(0.5)2  -  (0.0625)2 


=  0.19048 


from  Eq.(3 .46),  we  have 


k->  —  ^  —  0. 19048 

s2  =  b2,0  =  1 


(A. 8) 


(A.  9) 


The  new  polynomial  is 


H2(z)  =  b20  +  b21z-1  +  b22z 


y'2 
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I  .  r  r  r  «  .  »  .  r  *  *»*  -  •  .  •  »  ■  .  -  i 


H2(z)  -  1  +  0.476 19Z*1  +  0.19048Z*2 


for  m  =  2  and  i-  0,  we  have 
b->  a  -  frib') 


l  -  kV 


1  -(0.19048X0.19048) 

I  -  (0.19048/ 


(A.  10) 


(A.  11) 


and  finally 


bU  = 


b2.1  *  k2b2. 

1  k/ 


0.47619  -  (0.19048X0.47619) 

b,  ,  -  - - - ^ - —  «  0.40000 

1  -  (0.19048)“ 


From  Eq.(3.46),  we  have 


(A. 12) 


kl  “  bl,l  ®  0.40000 


S1  •  bl,0  ■  1 


(A.13) 


Therefore,  reflection  coefficients  and  s  coefficients  of  the  FIR  lattice  structure  are 
kj  *  0.40000 
k2  =  0.19048 
k3  -  0.0625 

sl  =  1 


s3  =  0.5 
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Example  A. 2  ( N  =  4:even) 

Obtaining  a  lattice  structure  for  the  linear  phase  FIR  transfer  function.  A  linear 
phase  FIR  transfer  function  is 

H3(z)  =*  0.154  4-  0.462Z'1  4-  0.462z‘2  4-  0.154z‘3  ,'A.U) 

Solution :  Eq.(A.14)  can  be  written  in  the  form  of  Eq.(3.35) 

h3(z)  =  a0  4-  ajz'1  4-  z'2(a1  +  a^'1)  (A.15) 

Comparing  Eqs.(A.14)  and  (A.15),  we  get  the  following  relationships. 


aQ  =  0.154 
aj  =  0.462 


(A.  16) 


To  determine  the  lattice  reflection  coefficients  of  the  corresponding  linear  phase  FIR 
lattice  structure  and  the  number  of  unit  delay,  we  use  Eq.(3.29)  which  is 


HN-l(z)  “  fM(z)  +  z'°  Gm(z> 


(A.  1 7) 


From  Eq.(A.14),  we  have  N  =  4(even),  the  order  of  the  polynomials  F^(z)  and  Gjyj(z), 
M  =  (N/2)-l  =  (4/2)-l  =  1,  and  the  number  of  unit  delays,  D  =  N/2  =  2.  Now, 
from  Eq.(3.37),  we  may  have  the  forward  prediction  error  transfer  function  as  follows: 


F j(z)  -  a0  4-  aj  z'1 

=  0.154  4-  0.462  z'1 


(A. 18) 


for  VI  =  1.  Eq.(A.  IS)  can  be  rewntten  as 


From  Eqs.(A.18)  and  (A.  19),  we  get  the  desired  reflection  coefficient  and  s  coefficient 
of  the  linear  phase  FIR  lattice  structure  as  follows: 

kl  =  bl,l  =  °-462 
S1  =  bl,0  "  °'154 

Example  A. 3  (N —  5\odd) 

Obtaining  a  lattice  structure  for  the  linear  phase  FIR  transfer  function.  A  linear 
phase  FIR  transfer  function  is 

H4(z)  =  0.15  -  0.45Z'1  +  0.36z‘2  -  0.45z'3  +  0.i5z'4  (A.20) 

Solution:  Eq.(A.20)  can  be  written  in  the  form  of  Eq.(3.30) 

H4(z)  =  a0  +  ajz*1  +  (l/2)a2z'2  +  z‘2  {(l/2)a2  +  ajZ*1  +  a0z*2}  (A.21) 

Comparing  Eqs.(A.20)  and  (A.21),  we  get  the  following  relationships. 


aQ  =  0.15 

aj  =  -0.45  (A. 22) 

a2  *  0.36 

To  determine  the  lattice  reflection  coefficients  of  the  corresponding  linear  phase  FIR 
lattice  structure  and  the  number  of  unit  delay,  we  use  Eq.(3.29)  which  is 

Hn_!(z)  -  Fm(z)  +  z'D  Gm(z)  (A. 23) 

From  Eq.(A.20),  we  have  N—  5(odd),  the  order  of  the  polynomials  F^(z)  and  G-yj(z), 
\1  =  IN-1V2  =  (5-0/2  =  2.  and  the  number  of  unit  delays.  D  =  (N-1V2  =  2.  Now. 
from  cq.(  3.32).  we  may  have  the  forward  prediction  error  transfer  function  as  fcilows: 


F2(z)  =  a0  +  aj  z'1  +  (1/2)  a2  z'2 
-  0.15  -0.45  z*1  +  0.18  z’2 


(A. 24) 


for  M  “  2,  Eq.(A.24)  can  be  rewritten  as 


.  it. 


*1  '4I  1 '4* 


g«  ,1.  il. 


F2(z)  =  b2>0  +  b2J  z'1  +  b2  2z'2 

From  Eqs.(A.24)  and  (A. 25),  we  have 

k2  *  b2  2  =  0.18 
$2  —  b2(Q  =  0.15 

Now,  we  need  to  generate  the  coefficients  for  Fj(z)  and  from  Eq.(3.46) 


(A.25) 


bl,0  ~ 


bU" 


s2b2.0  '  k2b; 


s2b2  1  *  k2b2 

— 1~T 

s2  -  &2 


(0. 15)(0. 15)  -  (0.18X0.18) 


(0.15)(-0.45)  -  (0.18)(-0.45)  , 

=  - - - - - =  -1.3ojo4 

0.152-  0.18" 


(A.26) 


(A. 27) 


I 


From  Eqs.(A.26)  and  (A.27),  we  have 
ki  =  bj  j  =  -1.36364 
S1  =  bl,0  =  1 

Therefore,  reflection  coefficients  and  s  coefficients  of  the  linear  phase  FIR  lattice 
structure  are 

kL  -  -1.36364 
k2  =  0.18 


s2  »  0.15 


rw 


APPENDIX  B 
COMPUTER  PROGRAMS 


rogrclin  1  '***7C*****’,f***  X)*f*XrCKXX*X*XXX7rX^XKXKyt*XK 

*  This  program  is  designed  for  the  system  identification  experiment  * 

*  which  is  shown  in  Section  (IV. D).  The  learning  curves  can  be  obtained* 

*  by  plotting  the  error,  E,  versus  the  update  iteration,  k.  * 

******** ****************** *********************************************** 


* 

* 

* 

* 

* 

* 

* 

* 

* 

★ 


* 

* 

* 

* 

* 

★ 

* 

* 

* 

* 

* 


Variable  Definition 

ISEED  :  seed  for  the  random  number  generation  (white  noise) 
AMU  :  adaptation  constant 

k  :  time  index 

M  :  order  of  the  FIR  transfer  function  or  total  number  of 

lattice  sections 
NB  :  number  of  iterations 

A(I)  :  coefficients  of  the  FIR  transfer  function 

B ( I )  :  reflection  coefficients  of  the  lattice 

F(I)  :  forward  prediction  error 

G( I )  :  backward  prediction  error 

GD(I)  :  delayed  backward  prediction  error 

SGL(I):  estimations  of  power 

W(K)  :  input  of  both  FIR  and  lattice 

YF (K)  •.  output  of  the  FIR  filter 

YL(K)  :  output  of  the  lattice  filter 

ER(K)  :  squared  error 

Variable  Declaration 


INTEGER  ISEED ,K,I,J,N,M,NB 

REAL  A(100) ,B(100),F(100),G(100),GD(100) ,YF(10000) ,YL(10000) 
REAL  X(100) .SGL(IOO) , ER( 10000 ) ,W( 10000) ,Y (10000) ,AMU,E 


Set  .^dictation  Constant  U  and  Number  of  Iterations  MB 


1 


* 


* 


FORMAT  (5X,  'AM’J'  ,5X,  '  NB  ’  ) 
READ (5,*)  AMU,NB 


Initialization 


..  ,,  ...  .»»  laU.)4lUUl,kCl^ldl,i 


B(K)=0 
X(K)=0 
F(K)=0 
G(K)=0 
GD(K)=0 
SGL(k)=l .0 
10  continue 

DO  15  K=l, 10000 
YF (K)=0 
YL(K)=0 
ER(K)=0 
W(K)=0 
Y(k)=0 

15  CONTINUE 
E=0. 

R=1 

ISEED=343169 

M=2 

A(  1 )  =  1.0 

A(2)  =  -0.89 
A(3 )  =  +0.25 

it 

k  Random  Number  Generation 

k  (mean  zero,  unit  variance,  white  sequence) 

k 

DO  20  N=1 ,NB 

CALL  LNORM(ISEED,RN, 1,1,0) 

W(N)  =  RN 
20  CONTINUE 

k 

k  FIR  Filter  Calculation 

* 

DO  30  K  =  1 ,NE 
X(1)=W(K) 

YF(K)=  A(1)*X(1) 

DO  40  1=1, M 

YF(K)=YF(K)+A( 1+1 )*X(I+1 ) 

40  CONTINUE 

DO  45  1=1, M 

X(M+2-:)=*X(M+I-I) 

45  CONTINUE 

30  CONTINUE 

A 

*  Lattice  Filter  Calculation 


DO  50  K  =1  ,NB 
F(1)=W(K) 


UVV  t.t. 


GU)«W<R) 

DO  60  I  ■  i  H 

P(I+1)»P(I)*B(I )*GC( I ) 

G(l*l /-GDC  ;*B( 

BO  CONTINUE 

/  ]r  \  >r  /  ' 

Calculating  tha  Error 
E  •  YF(K)  -  YL(K) 

Updating  tha  Raflaction  Coefficients 
CALL  LMS  f  9 . GD . E , AMU , SGL , M ) 

ER(K)-E**2 

DO  70  J-l.M 
GD(J)*G( J) 

70  CONTINUE 

IF  (K.EQ.R)  THEN 

WRITE ( h  ,  300  K .  9(1)  3 ' 2  ER'K. 

R*R+50 
END  IF 
Y(k)*E 

50  CONTINUE 

300  FORMAT (3X ,I6,5X,3(F10.7,4X)) 

Plotting  the  Learning  Curve 

CALL  PL0T(Y,N) 

STOP 

END 


SUBROUTINE  LMS ( B , GD , E , AMU , SGL , M ) 

REAL  B(100) ,GD(100) , SGL (100) ,E,AMU 
INTEGER 

go  ico 

SGL  ( 1 )  =0 . 3  'SGL  Vi,*':..';  CD  •  : ,  'GD  c ,  / 

CONTINUE 
DO  210  1=1 ,M 

B(I)=B{I)+(AMU/SGL(I))*E*GD(I) 

CONTINUE 

RETURN 


v.v’v.v 


SUBROUTINE  PLOT '  Y  .  N ) 

DIMENSION  Y < N ) ,X( 10000) 

ISTP-N/10 
DC  C  >1  N 
X( J  |*J 
CALL  TEK618 
CALL  PRTPLT (72,6) 

CALL  SHERPA(  PARKPARK ' . ' A ' , 3 ) 

CALL  PM7SCP ' l  1  ‘ 

CALL  RESET (  ALL  ) 

CALL  PAGE' 8.5, 11.0) 

CALL  HWROT I  AUTO ’ ; 

CALL  XINTAX 

CALL  AREA2D ( 5 . 0 . 2 . 8  J 

CALL  HE  I GHT '  C  .  1  C  ' 

CALL  COMPLY 

CALL  SHDCHRl 90. 0.1, 0.002,1) 

CALL  HEADIN'  'LEARNING  CURVES '  . 100, 2 .0 , 1 ) 

call  -c;ame  :ts?at::ns3  ::: 

CALL  YNAME  <  ERRORS  , 100 ) 

CALL  ME5SAG'  ADAPTIVE  LATTICE ( AMU* . 5 , FIG4 . 8 ) $ 1 ,100,3.0,-0.8) 

CALL  FRAME 

CALL  GRAF (0,1 STP ,N, *3. 0,1. 5, 3.0) 

CALL  THKCRV (0.02) 

CALL  CURVE (X,Y.N,0) 

CALL  ENDPL(O) 

CALL  DONEPL 

RETURN 

END 


& 


******p  rogram  2*****************,'r,lt,,[*******,lt*,’r,'r,'r*1*t***,*r***,,r**************** 

*  This  program  is  designed  for  the  system  identification  experiment  * 

*  using  the  LMS  algorithm  which  was  derived  in  Section  (V.B).  * 

************************************************************************* 


INTEGER  IS£ED,K, I , J,N,M,NB 

REAL  A(100) ,B(100) , F(1Q0) ,G{100) ,GD(100) ,x(100) ,YF(5000) ,YL(5000) 
REAL  ER(5000) ,W(50Q0) ,SGL(100) ,PH(100,100) ,PS(100,100) 

REAL  PSD(100 , 100) ,GR(100) ,Y(5000) ,AMU,E 

Set  Adaptation  Constant  Jl  and  Number  of  Iterations  NB 

WRITE (5,1) 

FORMAT ( 5X , 'AMU' ,5X, 'NB' ) 

READ ( 5 , *)  AMU,NB 

Initialization 

E=0. 

R=1 

DO  5  K=  1,5000 
W(K)=0 
YF(K)=0 
YL(K)*0 
ER(K)=0 
Y(K)=0 
CONTINUE 

DO  10  K=1 , 100 
A(K)=0 
B(K)=Q 
X(K)=0 
F(K)=0 
G(K)=0 
GD(K)=0 
SGL(K)=1 .0 
GR ( K ) =0 
CONTINUE 
DO  15  K-l , 100 
DO  15  L=1 , 100 
PH(K,L)=0 
PS(K,L)=0 
PSD(K,L)=0 
CONTINUE 
ISEED=343169 
M=2 


& 


|tt< 

I 

w 

v' 

ti1' 


B 


S 


A(l)  =  1.0 

A(2)  =  -0.89 
A(3)  =  +0.25 

Random  Number  Generation 
(mean  zero,  unit  variance,  white  sequence) 

DO  20  N=1 ,NB 

CALL  LNORM(ISEED,RN, 1,1,0) 

W(N)  =  RN 
20  CONTINUE 

FIR  Filter 

DO  30  K  =  1,NB 
X(1)=W(K) 

YF(K)=  A(1)*X(1) 

DO  40  1=1, M 

YF(K)=YF(K)+A(I+l)*X(I+l) 

40  CONTINUE 

DO  45  1=1, M 

X(M+2-I)=X(M+l-I) 

45  CONTINUE 

30  CONTINUE 

Lattice  Filter 

DO  50  K  =1  ,NB 
F(1)=W(K) 

G(1)=W(K) 

DO  60  I  =  1 ,M 

F(I+1)=F(I)+B(I) *GD ( I ) 

G(I+1)=GD(I)+B(I)*F(I) 

60  CONTINUE 

YL(K)=F(M+1) 

E  =  YF (K)  -  YL(K) 

( 

:  Updating  the  Reflection  Coefficients 

CALL  MLMS  (?H,?S,?SD,GR,F,G,GD,3,SGL,E,AMU,i1) 

ER(K)=E**2 
DO  70  J=1,M 
GD(J)=G(J) 

70  CONTINUE 

IF  (K.EQ.R)  THEN 

WRITE (6, 300)  K,  B(l) ,B(2) ,B(3) ,B(4) ,B(5) ,B(6) ,B(7) ,B(8) ,E 
R=R+30 


86 


END  IF 
Y(X)-£ 

c  WRITE(6 , 100)  K,W(K),YF(K),YL(K),E,ER(K) 
50  CONTINUE 

300  FORMAT ( 1X,I6,2X,9(F10.7,2X)) 

100  FORMAT (3X . 13 . 2X . 5(F1Q .7 , 2X} ) 

k 

k  Plotting  the  Learning  Curve 

k 

CALL  PLOT ( Y , N ) 

STOP 

END 


SUBROUTINE  MLMS ( PH , PS , PSD , GR , F , B , BD , R , SOL , EK , AMU , N ) 

DIMENSION  PH( 100 , 100 ) ,PS ( 100 , 100 ) , PSD (100 , 100 ) ,F(100) ,B<100) 
1 ,BD ( 100 ) , R( 100 ) , GR( 100 ) , SGL( 101 ) 

GR(N)*BD(N) 

DO ^ 200  1*2, N 
?Hf:,l>BD(N+l-H 
PS ( 1 , 1 )*F (N+l-I ) 

IT»I-1 

DO  10  K*1 , IT 

PH(I,K+1)*PH(I,K)+R(N+1-I+K)*PSD(I,K) 

PS( I , K+l )*PSD(I , K)+R(N+1-I+K) *PH(  I , K) 

DO  20  K*1 , IT 
PSD (I ,K)*PS(I ,K) 

CONTINUE 
DO  210  K*2  ,N 
GR(N+1-K)*PH(K , K) 

DO  211  K*1,N 

SGL(K)».90*SGL(K)+.10*GR(K)*GR(K)+1.0 
DO  220  1*1, N 

R(I)*R(I)+(AMU/SGL(I))*EK*GR(I) 

IF(R(I) .GE.1.0)  R(I )=0 , 

IF(R(I) .LE.-1.0)  R( I )=0 . 

CONTINUE 

RETURN 


SUBROUTINE  PLOT(Y,N) 
DIMENSION  Y(N) ,X(5000) 
ISTP*N/10 
DO  10  J=1,N 


%  \% 


^  »  j  *  •  **  ’ 


,  »*  4  aL  ^  M>  ,L  .L  .»  a  M 


X(J)-J 

CALL  TEK618 

CALL  PRTPLT (72,6) 

CALL  SHERPA(  PAAKPARK' .  A' ,3) 

CALL  PHYSOR ( 1 .  ,  1 .  ) 

CALL  RESET ( ' ALL ' ) 

CALL  PACE ( 3 . 5 , 1 i . 0 ) 

CALL  HWROT ( 'AUTO' ) 

CALL  XIHTAX 

CALL  AR£A2D( 5 . 0 . 2 .8) 

CALL  HEIGHT (0.10) 

CALL  COHPLX 

CALL  SHDCHR (90.0.1,0. C02 , 1 ) 

CALL  HEAD  IN  ( 'LEARNING  CURVES  MOO  .  2 .0 . 1 ) 

CALL  XNAME (' HERAT IONSS  ,100) 

CALL  YNAME ( ' ERRORS ' .100) 

CALL  HESSAG(  ADAPTIVE  LATTICE ( AMU- . 5 , FIG5 . S ) S 1 , 100  3. 0,-0. 8) 
CALL  FRAME 

CALL  GRAF(0 , ISTP ,N, -2. 5, 1.25.2.5) 

CALL  THKCRV(0.02 
CALL  CURVE  <  ?  N  O' 

CALL  ENDPL \ 0 ; 

CALL  DONE PL 

RETURN 

END 


A*****Pf Qgg-fjn  3********************************************************** 

*  This  program  is  designed  for  the  system  identification  experiment.  * 

*  The  LMS  algorithm  shown  in  Section  (V.C)  was  extended  to  the  linear  * 

*  phase  FIR  lattice  filter.  * 

**************** AAA A ************ ***************************************** 


INTEGER  I SEED, K, I ,J,N,M,NB,R,MA,ML 

REAL  A( 100) ,8(100) ,F< 100), G( 100) ,GD( 100), YF(3000; , YL(3000) 

REAL  X(100'.  H(  100 >  .YD (3000)  , ER( 3000)  ,W(3000) , YFF ( 3000)  , YLL(3000) 

REAL  SGL(IOO) , PH( 100 , 100 ) , PS( 100 , 100 ) , PSD( 100 , 100 ) 

REAL  GR( 100 ) , GRD( 100,100) , GRB( 100 ), Y( 3000) ,AMU,E , C 

Set  Adaptation  Constant  fl  and  Number  of  Iterations  NB 

WRITE (5,1) 

FORMAT (5X, 1  AMU ’ , 5X, ’ NB ’ ) 

READ ( 5 , *)  AMU ,NB 

Initialization 

E«0. 

R-l 

M«4 

ISEED-343169 
DO  10  K=1 , 100 
A(K)«0 
B(K)=0 
X(K)*0 
F(K)*0 
G(K)=0 
GD(K)=0 
H(K)=0 
SGL(K)=1.0 
GR(K)=0 
GRB(K)=0 
CONTINUE 
DC  11  X= 1,2000 
W(KJ=0 
YF(R)=0 
YL(K)=0 
YFF (K)=0 
YLL(K)=0 
YD(K)=0 
ER(K)*0 
Y(K)=0 


■t  i'  >1  ■!  . 
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CONTINUE 
DO  15  K*1 , 100 
DO  15  L=1 , 100 
PH(K,L)*0 
PS(K,L)=0 
PSD(K,L)=0 
GRD(K,L)=0 
CONTINUE 
C*.154 
A(l)*. 154/C 
A(2)a . 462/C 
C-.15 

A(l)*.15/C 

A(2)»-.45/C 

A(3)a.36/C 

Random  Number  Generation 
(mean  zero,  unit  variance,  white  sequence) 

DO  20  Nal,NB 

CALL  LN0RM(ISEED,RN, 1,1,0) 

W(N)  *  RN 
CONTINUE 

IF  (MOD(M,2) .EQ.0)  MA=M/2 
IF  (MOD(M,2) .NE.0)  MA=(M+l)/2 
IF  (MOD(M, 2) . EQ.0)  ML*MA 
IF  (MOD(M,2) .NE.0)  ML*MA-1 

Separation  of  Coefficients 

IF  (MOD(M, 2) .EQ.0)  THEN 
DO  21  1*1, MA 
H(I)*A(I) 

H(MA+2+I)*A(MA+l-I) 

CONTINUE 

H(MA+l)=A(MA+l)/2 
H(MA+2)=A(MA+l)/2 
END  r? 

:?  !M0D(M.Z) .ne.o)  then 

DO  22  :  =  L,:iA 

H(I)=A(I) 

H(MA+I)*A(MA+1-I) 

CONTINUE 
END  IF 

LINEAR  PHASE  FIR  FILTER 


K 


% 


*.  *'  - 


%  %  v-  *'■  ■. 


■ 


DO  30  K  *  1,NB 
X(1)=W(K) 

YF(K)=H(1)*X(1) 

IF  (MOD(M,2).EQ.O)  THEN 
DO  40  1=1,  MA 

YF(K)=YF(K)+H(I+1)*X(I+1) 

40  CONTINUE 

YFF(K)=YF(K) 

DO  41  1=1 ,MA+1 

YFF (K)=YFF(K)+F(MA+1+I)*X (MA+I) 

41  CONTINUE 
END  IF 

IF  (MOD(K,2).NE.O)  THEN 
DO  42  1*1 ,MA-1 

YF(K)=YF(K)*H(I+1)*X{I+1) 

42  CONTINUE 

YFF (K)=YF (K) 

DO  43  1=1, MA 

YFF ( K ) =YFF ( K ) +H ( MA+I ) *X ( MA+I ) 

43  CONTINUE 
END  IF 

DO  45  1*1, M 

X(M+2-I)=X(M+l-I) 

45  CONTINUE 

30  CONTINUE 

r 

'  LATTICE  FILTER 

r 

DO  50  K  =1 ,NB 
F(1)=W(K) 

G(1)=W(K) 

DO  60  I  =  1 ,ML 

F(I+1 )=F( I )+B ( I )*GD( I ) 

G(I+1)=GD(I)+B(I)*F(I) 

60  CONTINUE 

YL(K)=F (ML+1 ) 

YD(K)=G(t!L+l) 

YLL ( K ) =YL ( K ) +YD  < K- MA ) 

E  =  YFF (IV,  -  YLL ,  K ) 

*  Updating  the  Reflection  Coefficients 

« 

CALL  MLMS  ( P , Q , QD , GRB , GRD , F , G , GD , B , SGL , E , AMU , ML , MA ) 
ER(K)=E**2 
DO  70  J=1,ML 
GD{J)=G(J) 

70  CONTINUE 


91 
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Y(K)=E 

WRITE(6 ,300)  K,Y(K) 

50  CONTINUE 

it 

k  Plotting  the  Learning  Curve 

lr 

CALL  PL0T(Y,N) 

300  F0RMAT(3X,I6,3X,F10.5) 

STOP 

END 


SUBROUTINE  MLMS ( PH , PS , PSD , GRB , GRD , F , B , BD , R , SGL , EK , AMU , N , MA ) 
DIMENSION  PH(100, 100) , PS (100 , 100) , PSD (100 ,100) , GRD (100 ,100) 
1 , GRB (100) ,F(100) ,B(100) ,BD(100) ,R(100) ,GR(100),SGL(101) 
GR(N)=BD(N) 

GRB(N)=F(N) 

DO  200  1=2, N 
PH(I,1)=BD(N+1-I) 

?S(:,1)=F(N+1-I) 

IT=I-1 

DO  110  K=1 , IT 

PH(I.K+1)=PH(I,K)+R(N+1-I+K)*PSD(I,K) 

PS(I,K+1)=PSD(I,K)+R(N+1-I+K)*PH(I,X) 

DO  120  K=1 , IT 
PSD (I ,K)=PS(I ,K) 

CONTINUE 

DO  210  K=2,N 

GR(N+1-K)=PH(K,K) 

GRB(N+1-K)=PS(K,K) 

DO  220  K=1,N 
DO  220  L=1,MA 

GRD(K,MA+2-L)=GRD(K,MA+l-L) 

DO  230  K=1,N 
GRD(K,1)=GRB(K) 

DO  240  K=1,N 

GRB ( K ) =GRD < K , MA+1 ) 

DO  Z50  K=i,N 

gr(x)=gr(,k;+grb(x) 

DO  260  K=1,N 

SGL(K)=.90*SGL(K)+.10*GR(K)*GR(K)+1 .0 
DO  270  1=1, N 

R ( I ) =R ( I ) + ( AMU/ SGL ( I ) ) *EK*GR ( I ) 

IF(R(I) .GE.1.0)  R(I)=0. 

IF(R(I) .LE.-1.0)  R(I)=0 . 

CONTINUE 


i 


K£ 


B 


RETURN 

END 

* 

* 

★ 

SUBROUTINE  PLOT(Y,N) 

DIMENSION  Y(N) ,X(3000) 

ISTP=N/10 
DO  10  J=1,N 
10  X(J)*J 
C  CALL  TEK618 

c  CALL  PRTPLT(72 , 6 ) 

CALL  SHERPA( ' PARK? ARK 1 , 1 A ' , 3 ) 

CALL  PHYSOR(l. .1.) 

CALL  RESET ( 'ALL' ) 

CALL  PAG£(8.5.11.) 

CALL  HWROT( ' AUTO' ) 

CALL  XINTAX 

CALL  AREA2D(5.0,2.8) 

CALL  HEIGHT(O.IO) 

CALL  COMP LX 

CALL  SHDCHR(90. 0.1, 0.002,1) 

CALL  HEAD IN ( 'LEARNING  CURVES ', 100 , 2 . 0 , 1 ) 

CALL  XNAME ( ' ITERATIONSS ' , 100 ) 

CALL  YNAME ( ' ERRORS ',100 ) 

CALL  MESSAG( 'ADAPTIVE  G-M  LATTICE  (F5 . 10  ,  AMU*-0 .1)$'  ,100,3.0,-0.8 
CALL  FRAME 

CALL  GRAF(0 , ISTP ,N , -8. 0,4. 0,8.0) 

C  CALL  THKCRV(0 . 02 ) 

CALL  CURVE (X,Y,N,0) 

CALL  ENDPL(O) 

CALL  DONEPL 

RETURN 

END 


*  This  program  is  designed  for  the  estimation  of  spectral  lines  ; 

*  white  noise.  The  input  process  xvk)  consists  of  a  signal  in  nc; 

*  and  a  signal  may  be  a  single,  or  multiple  sinusoids.  The  algci ; 

*  was  derived  in  Section  (V.D). 


INTEGER  ISEED.K.I . J.N.M.NB.R.MA  ML  D 

REAL  A(100)  ,B(I00dF'100)  ,G  100)  GD  i  190  '  YL  < 5000  ; 

REAL  Y(  100)  .  YD i  5000  )  .  W^SOOO  ,  .  YLLdOOO  >.  YD<  5000  '  IN?.  SJ 
REAL  SGLc 101 ) . PH; IOC . IOC ; , PS  100, 101  •  PSD  1O0.::: 

REAL  GRD(  100.  IOC'  GRBUOO)  .GR  100) 

REAL  RE  (1 00  MM.  100)  .AJ(iOO)  SD  .  SNR  AVG  AMP  AMU  E 


Set  Adaptation  Constant  ^  and  Number  of  Iterations  NE 


WRITE (5,1) 

FORMAT ( 5X ,  AMU  1  . 5X  1 NB  ) 
READ; 5  *)  AMU  NE 


Initialization 


ISPEC*1000 
SNR-30 . 

SD=10**(-(SNR/20)) 
AMP=*SQRT(2.  ) 

AVG*0 . 

E*0 . 

R=1 

PI=3. 141592634 
M=3 

MP1=M+1 

ISEED=343169 

IF  (MOD(M,2) .EQ.O)  THEN 

MA=M/2 

ML=MA 


ML=MA-1 
END  IF 

DO  10  K=1 , 100 
A(K)=0 
B(K)=0 
F(K)=0 
G(K)=0 


WiWWFiW 


GD(K)=0 

SGL(K)=1.0 

GR(K)=0 

GRB(K)=0 

RE(K)=0 

IM(K)*0 

AJt,K)=0 

Y(K)=0 

10  CONTINUE 

DO  11  K=1 , 5000 
W(K)=0 
YL ( K ) =0 
YLL(K)=0 
YD(K)=0 
ER(K)=0 
INP (K)=0 

11  CONTINUE 

DO  15  K=1 , 100 
DO  15  L=1 , 100 
PH(K,L)=0 
PS(K,L)=0 
PSD(K,L)=0 
GRD(K,L)-0 

15  CONTINUE 

A 

*  Random  Number  Generation 

*  (mean  zero,  unit  variance,  white  sequence) 

* 

DO  20  N=1,NB 

CALL  LNORM(ISEED,RN, 1,1,0) 

W(N)  =  SD*RN+AVG 

20  CONTINUE 

* 

*  LATTICE  FILTER 

A 

DO  50  K  =1  ,NB 
AK=K-1 

c  INP(K)=AMP*COS(2*?I*.15*AK)+W(K) 

-  :np(x;=amp*< cds^pi* 

ZNP  ( K )  SAMP’*  ( COS  ( 2*?I  *  -  2  5  *  AK  ‘CCS. 

*)+COS(2*PI*.35*AK))+W(X) 

F(1)=INP(K) 

G(1)=INP(K) 

DO  60  I  =  1,ML 

F(I+1)=F(I)+B(I)*GD(I> 
G(I+1)-GD(I )+B( I )*F ( I ) 

60  CONTINUE 


AD-A184  97* 
UNCLASSIFIED 


AN  ADAPTIVE  LATTICE  ALGORITHM  FOR  SPECTRAL  LINE 
ESTIHATION(U)  NAVAL  POSTGRADUATE  SCHOOL  HONTEREV  CA 
I  K  PARK  JUN  87 

F/Q  12/1 


J 


,*« 

% 


YL(K)»F(ML+1) 

YD(K)*G(ML+1) 

YLL(K)*YL(K)+YD(K-MA) 

E  *  YLL(K) 

Updating  the  Reflection  Coefficients 

CALL  MLMS  (PH,PS,PSD,GR,GRB,GRD,F,G,GD,B,SGL,E,AMU,ML,MA) 
ER(K)=E**2 
DO  70  J=1 ,ML 
GD ( J ) =G { J ) 

CONTINUE 

IF  (R.NE.ISPEC)  GO  TO  50 
IF  (K.EQ.NB)  THEN 
WRITE (6 ,300)  K, INP(K) ,E 

WRITE (6, 301)  K,B(1),B(2),B(3),B(4),B(5),B(6),B(7),B(8),B(9) 
R=R+100 

Determine  the  FIR  Coefficients  from  the  Lattice  Reflection 
Coefficients 

CALL  STEPUP  (A,B ,ML) 

IF  (MOD(M, 2) .EQ.0)  THEN 
DO  60  1*1 ,M/2 
A(M+2-I)=A(I) 

A(M/2+l )-2*A(M/2+l ) 

ELSE 

DO  81  I=l,(M+l)/2 
A(M+2-I)=A(I) 

END  IF 

WRITE  (6,600)  (I ,A(I) ,1=1 ,MP1) 

Calculate  the  Power  Spectrum 

CALL  SPEC  (RE,IM,A,AJ,M,PI) 

D=100 

WRITE  (6,601)  ( .  Q05*Z,AJ (Z) ,Z=1 ,D) 

Plotting  the  Spectrum 

CALL  PL0T(AJ,D) 

ISPEC=ISPEC+1000 
END  IF 

CONTINUE 

WRITE  (6,300)  (YLL(K) ,K=1 ,NB) 

FORMAT(10(F10 .7 , IX) ) 

F0RMAT(1X, 15 , IX, 10(F10 . 7 , IX) ) 


IT" 

a 

w 

$ 

a 

u 

! 


% 


’I 

$ 

£ 

i.. 

$ 

$ 

k. 

£ 


I 

I 


Vi 

R 

tfi 

••S' 

k 

ha 

rk» 

$ 

$ 

i 


& 

as 

& 

600  FORMAT  (1X,I5,5X,F15.7) 

601  FORMAT  (1X,F5.3,5X,F15.7) 

STOP 

END 

* 

* 

* 

SUBROUTINE  MLMS (PH , PS , PSD , GR , GRB , GRD , F , B , BD , R , SGL , EK , AMU , N , MA) 
DIMENSION  PH(100,100) ,PS(100,100) ,PSD(100 , 100) ,GR(100) ,GRB(100) 
1 , GRD (100 , 100) ,F(100) ,B(100) ,BD(100) ,R(100) ,SGL(101) 

GR(N)=BD(N) 

GRB(N)=F(N) 

DO  200  1=2,  N 
PH(I , 1)=BD(N+1-I) 

PS(I,1)=F(N+1-I) 

IT=I-1 

DO  110  K=1 , IT 

PH (I ,K+1)=PH(I ,K)+R(N+1-I+K)*PSD(I (K) 

PS (I ,K+1)=PSD(I ,K)+R(N+1-I+K)*PH(I ,K) 

110  CONTINUE 

DO  120  K=1 , IT 
120  PSD(I ,K)=PS (I ,K) 

200  CONTINUE 

DO  210  K=2 ,N 
GR(N+1-K)=PH(K,K) 

GRB(N+1-K)=PS<K,K) 

210  CONTINUE 

DO  220  K=1  ,N 
DO  220  L=1,MA 

220  GRD(K,MA+2-L)=GRD(K,MA+l-L) 

DO  230  K=1 ,N 
230  GRD(K, 1)=GRB(K) 

DO  240  K=1 ,N 
240  GRB ( K ) =GRD ( K , MA+ 1 ) 

DO  250  K=1 ,N 
250  GR(K)=GR(K)+GRB(K) 

DO  260  K=1 ,N 

260  SGL(K)=.90*SGL(K)+.10*GR(K)*GR(K)+1.3 

DO  270  1  =  1. M 

R(I)=R(I)-(AMU/SGL(I))*EK*GR(I; 

IF(R(I) .GE.1.0)  R(I )=0 . 

IF(R(I) .LE.-1.0)  R(I)=0. 

270  CONTINUE 
RETURN 
END 


SUBROUTINE  STEPUP  (A,B,ML) 
DIMENSION  A(IOO) ,C(100) ,B(100) 
A(l)*l. 

A(2)=B(1) 

DO  30  MINC=2,ML 
DO  10  J=1,MINC 
JB*MINC-J+1 
C(J)=A(JB) 

DO  20  IP=2 ,MINC 
A(IP)=A(IP)+B(MINC)*C(IP-1) 
A(MINC+1 )=B(MINC) 

CONTINUE 

RETURN 

END 


SUBROUTINE  SPEC(RE,IM,A,AJ,M,PI) 

REAL  RE(100) , IM(100) ,A(100) ,AJ<100) ,Y(100) 

RE  ( 1 ) =A ( 1 ) 

IM(1)=0, 

DO  91  J=l,100 
DO  92  1=1, M 

RE(I+1)=RE(I)+A<I+1)*C0S(2*I*PI*.5*J/100) 

IM(I+1)«IM(I)+A(I+1)*SIN(2*I*PI*.5*J/100) 

CONTINUE 

AJ(J)=-10.*ALOG10(RE(M+1)*RE(M+1)+IM(M+1)*IM(M+1)) 

CONTINUE 

TEMP=AJ ( 1 ) 

DO  93  L=2 , 100 

IF  (AJ(L).GT.TEMP)  TEMP=AJ(L) 

CONTINUE 

DO  94  L=1 , 100 

AJ(L)=AJ(L)-TEMP 

CONTINUE 

RETURN 


SUBROUTINE  PLOT(Y,N) 

DIMENSION  Y(N) ,X(100) 

ISTP=N/10 

DO  10  J=1,N 

X(J)=J 

DF=. 5/N 


DO  10  K*2,N 

X(K)»X(K-1)+DF 

XMIN=X(1) 

XMAX*X(N) 

XSTP»10*DF 
IYMIN=Y ( 1 ) 

IYMAX=Y(1) 

DO  20  K=2,N 

IF(Y(K) .GT. IYMAX)  IYMAX=Y(K) 

IF(Y(K) .LT.IYMIN)  IYMIN=Y(K) 

CONTINUE 

I YSTP= ( IYMAX - I YMIN ) / 5 

CALL  TEK618 

CALL  PRTPLT(72,6) 

CALL  SHERPA ( 1 PARKPARK 1 , 1 A 1 , 3 ) 

CALL  RESET ('ALL') 

CALL  PAGE (8. 5, 11.0) 

CALL  HWROT ( 1  AUTO 1 ) 

CALL  XINTAX 

CALL  AREA2D ( 5 . 0 , 2 . 8 ) 

CALL  HEIGHT(O.IO) 

CALL  COMPLX 

CALL  SHDCHR(90 .0 , 1 , 0.002 , 1) 

CALL  HEADINC FREQUENCY  SPECTRUMS ',100 , 2 . 0 , 1 ) 

CALL  XNAME  (  '  FREQUENCY$  MOO ) 

CALL  YNAME( 'MAGNITUDE (DB)$ 1 , 100) 

CALL  MESSAGC FIGURE  5.21  (M-8  ,SNR-30DB)$ ',100 ,3 .0 , -0 .8) 

CALL  MESSAG('MODEL  ORDER  SELECTION^  SINUSOIDS)$',100,3.0,-0.8) 
CALL  THKFRM(0.03) 

CALL  FRAME 

CALL  GRAF(XMIN,XSTP ,XMAX, IYMIN, IYSTP , IYMAX) 

CALL  THKCRV(0.02) 

CALL  CURVE(X,Y,N,0) 

CALL  ENDPL(O) 

CALL  DONEPL 

RETURN 

END 
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